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I Iniroduction 

This research was undertaken in an effort to aid the Forestry Application 
Project on Timber Resources. Mission M230 of the C~130 aircraft was flown over 
the Sam Houston National Forest on March 21, 1973 at 10,000 feet altitude. The 
Bendix 24 channel multispectral scanner collected the data. Four forest scenes of 
this data set were selected for study. They were edits 3, 6, 9, and 14. The cate" 
gories of timber classes and subclasses are shown in Table I.l . 

The application oriented research was to apply and document the capability 
of existing textural and spatial automatic processing techniques at the University of 
Kansas to classify the MSS imagery into specified timber categories. The ground 
truth for the study was supplied by the Forestry Applications Project. 

Over a hundred classification experiments were performed on this data 
using feature selected from the spectral bands and a textural transform band. The 
textural transform band is an image whose resolution cells have grey tone intensities 
which indicate one parameter of local neighborhood texture. The textural transform 
concept is discussed in Section III. The classification was done by equal interval 
quantizing the images to 32 levels and using a non~parametric table !ook"up rule 
discussed in Section II. The various spatial pre" and post"processing options are 
discussed in Sections IV and V. Sections VI through IX discuss the results using 
only spectra! features. Sections X through XIII discuss the combined spectral tex- 
tural results. 

The results indicate that 

(1) spatial post-processing a classified image can cut the classification 
error to'1/2 or 1/3 of its Initial value. 

(2) spatial post-processing the classified image using combined Spectral 
and textural features produces a resulting image with less error than 
post-processing a classified image using only spectral features. 

(3) classification without spatial post processing using the combined 
spectral textural features tends to produce about the same error rate 
as a classification without spatial post processing using only spectral 
features. 
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TABLE I.l THE TYPE (CLASSES) AND CONDITION CLASSES (SUBCLASSES) 

OF FOREST FEATURES OF INTEREST IN SAM HOUSTON NATIONAL FOREST OF TEXAS 


Type No. 

Type (Class) 

Subclass No. 

Condition Class (Subclass) 

1 

Shortleaf pine 

1.1 

Plantation “ 3 years old 



1.2 

Poletimber - immature 



1.3 

Sawtimber - immature 



1.4 

Sawtimber ~ mature 

2 

Loblolly pine 

2.1 

Plantation - 1 year old 



2.2 

Plantation - 3 years old 



2.3 

Seedling and Sapling - 




adequately stocked 

1 

. 

2.4 

Poletimber - immature 

1 


2.5 

Sawtimber - immature 

i 


2.6 

Sawtimber “ mature 

3 

\ 

Laurel oak - 

3.1 

Sawtimber - immature 


willow oak 



4 

Sweetgum ” nuttal 

; 4. 1 

•j 

Sawtimber - low quality 


oak ” willow oak 

j 4.2 

Sawtimber - immature 



1 4.3 

\ 

Sawtimber - mature 

5 

Post oak “ black oak 

5.1 

Sawtimber - immature 

6 

1 Loblolly pine “ 

6.1 

Sawtimber - immature 

1 

1 hardwoods 

1 


1 ^ 

1 Cut-over land 

7.1 

Site prepared and 


I windrowed 

I 7.2 Not site prepared 

t 

( 





These results mean that regardless of how the image is classified, spatial post”” 
processing should be used to reduce the error rate. Furthermore, the best 
post processing results can be obtained if textural features are used; but, if no 
spatial post“processing is going to be utilized, spectral bands only will give about 
the same results as the combined spectral textural bands. 

These conclusions are based on classification into all timber subclasses using 
large training sets averaging more than 25,000 points per image. Because the 
training sets were orders of magnitude larger than the number of categories times the 
number of features, the statistics must be considered as large sample statistics and 
we used, justifiably, the training data as the test data. 

Tables 1.2 and 1.3 summarize the basis of our conclusions. The results of 
each experiment can be summarized in three ways: by average error, by average 
misidentification error, and by average false identification error. The average 
error is defined as the total number of incorrect category assignments divided by 
the total number of assignments. The average misidentification error is defined as 
the equally weighted average over all categories of the number of times the category 
is incorrectly assigned divided by the total number of times the category occurs in the 
ground truth. The average false identification error is defined as the equally weighted 
average over all categories of the number of times an incorrect assignment is made to 
the category divided by the total number of times an assignment is made to the category 

When the ground truth has each category occurring with equal frequency, the 
average misidentification error will equal the average error. When the number of 
assignments to each category is the same, the average false identification error will 
equal the average error. If the prior probability for a category is high and the cate- 
gory has a high misidentification error, then all other things being equal, the average 
error will be higher than the average misidentification error. If the prior probability 
for a category is low, and the category has high misidentification error, then all 
other things being equal, the average error will be lower than the average misiden- 
tification error. 

From Tables 1.2 and 1.3 it is readily apparent that both the use of textural 
features and spatial post processing tends to increase and equalize the average mis- 
identification error and false identification error while cutting the average error 
to less than half its initial value. 
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1. 1 Contingency Tables of Classification Results 

All results ore reported with a complete contingency table. The contin" 
gency tables ore all organized in the same manner. The title for the contingency 
table tells which images are being compared. The first nine character file name 
is the name of the ground truth image file. The number following it is the sym- 
bolic band number used from that multi-image file. The second nine character 
file name is the name of the classified image file. The number following it is the 
number of the symbolic band used from that multi-image file. The row label 
UNKWN means unknown true category identification. The column label R DEC 
means reserved decision. 

The contingency tables have a column labeled ERR. This column designates 
the number of the resolution cells in each category misidentified . The next column 
is labeled % ERR and it designates the percent of misidentification error. The con- 
tingency tables have a row labeled ERR. This row designates the number of resolu- 
tion cells in each category fahely identified. The next row is labeled % ERR and 
it designates the percent of false identification error. The label % SD stands for 
the percent standard deviation of the error estimates. The entry whose row is 
labeled TOTAL and whose column is weighted % ERR is the equally weighted 
average of the misidentification error percentages. The entry whose column is 
labeled total and whose row is weighted % ERR is the equally weighted average of 
the false identification error percentages. 
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Table 1.2 summarizes the error rates obtained from the spectral versus the spectral- 
textural classification using 3 band pairs and no spatial post processing. 


Spectral Spectral-Texture 


t 

1 

I 

1 

Average i 
Error 

Average 
Misidenti- 
! fication 
! Error 

Average 
j False 
1 identifi- 
s cation 
' Error 

Average 
j Error 

Average 

Misidenti" 

fication 

Error 

Average 

False 

Identifi" 

cation 

Error 

Edit 6 i 

t 

9.3% 

■ 34% 

4— 

‘ 33% 

i 

1 6.8% 

38% 

1 

37% 

Edit 9 i 

! 19% 

i 25% 

; 32% 

j 

’ 15% 

27% 

33% 

Edit 14 

12% 

; 32% 

I 31% 

j 

‘ texture b 

and not selec 

ted by feature selector 

Edit 3 

24% 

1 35% 

1 40% 

i 12% 
1 

j 40P/o 

44% 


Table 1.3 summarizes the error rates obtained from the spectral versus the spectral 
textural classification using 3 band pairs and spatial post processing . 
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Table Look~Up Decision Rule 

Brooner, Haralick and Dinslein (1971) used a table look-up approach on 
high altitude multiband photography flown over Imperial Valley, California to 
determine crop types. Their approach to the storage problem was to perform an 
e qual probability quantizing from the original 64 digitized grey levels to ten 
quantized levels for each of the three bands; green, red, and near infrared. Then 
after the conditional probabilities were empirically estimated, they used a Bayes 

3 

rule to assign a category to each of the 10 possible quantized vectors in the 3- 
dimensional measurement space. Those vectors which occurred too few times in 
the training set for any category were deferred assignment. 

The rather direct approach employed by Brooner et a! . has the disadvantage 
of requiring a rather small number of quantized levels. Furthermore, it cannot be 
used with measurement vectors of dimension greater than four; for if the number of 
quantized levels Is about 10, then the curse of dimensionality forces the number of 
possible quantized vectors to an unreasonably large size. Recognizing the grey 
> level precision restriction forced by the quantizing coarsening effect, Eppier, 
Heimke, and Evans (1971) suggest a way to maintain greioter quantizing precision by 
defining a quantization rule for each category - measurement dimension as follows: 

(1) fix a category and a measurement dimension component; 

(2) determine the set of all measurement patterns which would be 
assigned by the decision rule to the fixed category; 

(3) examine all the measurement patterns in this set and determine 
the minimum and maximum grey levels for the fixed measurement 
component; 

(4) construct the quantizing rule for the fixed category and measurement 
dimension pair by dividing the range between the minimum and maximum 
grey levels for the category into equal spaced quantizing intervals. 

This multiple quantizing rule in effect determines for each category d 
rectangular parallelepiped in measurement space which contains all the measurement 
patterns assigned to it. Then as shown in Figure 11*1, the equal interval quantizing 
lays a grid over the rectangular parallelepiped. Notice how for a fixed number of 
quantizing levels, the use of multiple quantizing rules in each band allows greater 
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grey level quantizing precision compared to the single quantization rule for each 
band . 

A binary table for each category can be constructed by associating each 
entry of the table with one corresponding cel! in the gridded rectangular parallel” 
epiped. An entry is a binary 1 if the decision rule assigns a majority of the measure” 
ment patterns in the corresponding cell to the specified category; otherwise, the 
entry is assigned to be a binary 0. 

The binary tables are used in the implementation of the multiple quanti” 
zation rule table look~up in the following way. Order the categories in some 
meaningful manner such as by prior probability. Quantize the multispectral 
measurement pattern using the quantization rule for category c-j . Use the quan” 
tized pattern as an address to look up the entry in the binary table for category c^ 
to determine whether or not the pre-stored decision rule would assign the pattern 
to category c^ . If the decision rule makes the assignment to category c^ the entry 
would be a binary 1 and, all is finished. If the decision rule does not make the 
assignment to category c-j, the entry would be a binary 0 and the process would 
repeat in a similar manner with the quantization rule and table for the next category. 

One advantage to this form of the table look”up decision rule is the flexibility 
to use different subsets of bands for each category look”up table and thereby take 
full advantage of the feature selecting capability to define an optimal subset of bands 
to discriminate one category from all the others. A disadvantage to this form of the 
table look-up decision rule Is the large amount of computational work required to 
determine the rectangular parallelepipeds for each category and the still large amount 
of memory storage required (about 5,000 8 bit bytes per category). 

Eppler (1974) discusses a modification of the table look“up rule which enables 
memory storage to be reduced by five times and decision rule assignment time to be 
decreased by 2 times. Instead of pre -storing in tables a quantized measurement space 
image of the decision rule, he suggests a systematic way of storing In tables the 
boundaries or end-points for each region in measurement space satisfying a regularity 
condition and having all its measurement patterns assigned to the same category. 

Let D ^ D| X X. . .x be measurement space. A subset 

RC D| x D 2 X. . .X D|^ is a regular region if and only if there exists constants 
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and H| and functions L2/ . . . , H2, . . . , Hj^ 

((L^: X D2 X. . .X (-%“); H^: x D2 x. . .x -> (-“,“>)) 

such that 

R = |(^1 1 • • • r 


' ^2' ' ’ ■ ^ ' ^ 2 ' ' ^N"l 

From the definition of a regular region, it is easy to see how the table look-up 
by boundaries decision rule can be implemented. Let d = (d^, . . ., dj^) be the measure- 
ment pattern to be assigned a category. To determine if d lies within a regular region R 
associated with category c we look up the numbers L-j and H| and test to see if d^ lies 
between and . If so, we look up the numbers L 2 (di) and H 2 (di) and so on. If 
all the tests are satisfied, the decision rule can assign measurement pattern d to category 
c. If one of the tests fails, tests for the regular region corresponding to the next category 
can be made . 

The memory reduction in this kind of table look-up rule is achieved by only 
storing boundary or end-points of decision regions and the speed-up is achieved by 
having one-dimensional tables whose addresses are easier to compute than the three 
or four-dimensional tables required by the initial table look-up decision rule. However, 
the price paid for by these advantages is the regularity condition imposed on the decision 
regions for each category. This regularity condition is stronger than set connectedness 
but weaker than set convexity. 

Another approach to the table look-up rule can be based on Ashby's (1964) 
technique of constraint analysis. Ashby suggests representing in an approximate way 
subsets of Cartesian product sets by their projections on various smaller dimensional 
spaces. Using this idea for two-dimensional spaces we can formulate the following 
kind of table look-up rule. 

Let D = X D 2 X. . .X D|^ be measurement space, C be the set of categories, 
and J C {1, 2,..., Nl x {1,2,..., N)be an index set for the selected two-dimensional 


)} 
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spaces. Let the probability threshold a be given. Let (i,j)^j; for each X2)«^Dj x D 

define the set S..(x.j, x^) of categories having the highest conditional probabilities 
ij I z 

given (x^ , by 


(c) > a 


S..(x-|, Xr,)- {ccCIp (c) ^*a..)}, where a., is the largest number 

IJ I Z X-j,X^ Ij Ij 

which satisfies 

E p 

ceS..(x,,x.) 

S..(xi, X2) is the set of likely categories given that components i and j‘ of the 

measurement pattern take the values (x^, X2). 

The sets S.^., (i,j)eJ, can be represented in the computer by tables. In 

the table S.j the (x-j, ^2/^ entry contains the set of all categories of 

sufficiently high conditional probabilities given the marginal measurements 
(xp x^) from measurement components i and ], respectively. This set of categories 
is easily represented by a one word table entry; a set containing categories Cy, 
c^, and c^2/ f°*" example, would be represented by a word having bits 1 , 7 , 9 , and 
12 on and all other bits off. 

The decision region R(c) containing the set of all measurement patterns to 
be assigned to category c can be defined from the S.j sets by 

R(c) = |(d^, d2, . . ., dj^)eD^ X D2 X. . .X D|^! {c} = n S..(d.,d.)| 


This kind of a table look~up rule can be implemented by using successive pairs of 
components (defined by the index set J) of the (quantized) measurement patterns 
as addresses in the just mentioned two-dimensional tables. The set intersection 
required by the definition of the decision reg ion R(c) is implemented by taking the 
Boolean AND of the words obtained from the table look-ups for the measurement to be 
assigned a category. Note that this Boolean operation makes full use of the natural 
parallel compute capability the computer has on bits of a word. If the k bit is the 
only bit which remains on in the resulting word, then the measurement pattern is 
assigned to category Cj^. If there is more than one bit on or no bits are on, then 
the measurement pattern is deferred its assignment (reserved decision). 
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Thus we see that this form of a table look-up rule utilizes a set of "loose" 
Bayes rules in the lower dimensional projection spaces and intersects the resulting 
multiple category assignment sets to obtain a category assignment for the measure- 
ment pattern in the full measurement space . 

Because of the natural effect which the category prior probabilities have 
on the category assignments produced by a Bayes rule it is possible for a measure- 
ment pattern to be the most probable pattern for one category yet be assigned by 
the Bayes rule to another category having much higher prior probability. This 
effect will be pronounced in the table look-up rule just described because the 
elimination of such a category assignment from the set of possible categories by 
one table look“up will completely eliminate it from consideration because of the 
Boolean AND or set intersection operation. However, by using an appropriate 
combination of maximum likelihood and Bayes rules, something can be done about 
this. 


For any pair (i,j) of measurement components, fixed category c; and pro- 
bability threshold &, we can construct the set of Tjj(c) having the most probable 
pairs of measurement values from component i and j arising from category c. The 
set '5 defined by 


Tjj(c) — ^ ^j 1 ^ 2 )^ 


where 3 jj(c) is the largest number satisfying 


^ 3 

(x^,X2)eT..(c) 

Tables which can be addressed by (quantized) measurement components 
can be constructed by combining the 5 .. and T.. sets.. Define Q.j(x.|, X2) by 


Qjj(x-j , ^2) IccC j (x-j , X 2 )€T.j(c)|uS..(x| , X2) 


The set Q..(x^, X2) contains all the categories wfiose respective conditional pro- 
babilities given measurement values (xp X2) of components i and j are sufficiently 
high (a Bayes rule criteria) as well as all those categories whose most probable 
measurement values for components i and j respectively are (xy , X2) (a maximum 
likelihood criteria). A decision region R(c) containing all the (quantized) measure- 
ment patterns can then be defined as before using the Q.. sets: 
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A majorfJ-y vote version of this kind of table look“up rule can be defined 
by assigning a measurement to the category most frequently selected in the lower 
dimensional spaces. 

R(c) — |(d I /^2' • • * / 1 ^ 1^2 

^ |(i,j)fj|ceQ.^(d.,d.)| ^ **^|(i,j)eJlceQ..(d.d.) 
for every ceC ~ (c }| 

Classification results were run with 3 = .07a and a chosen to minimize the 
number of reserved decisions. Figure II. 2 illustrates a graph of the number of 
reserved decisions versus probability threshold a. 



igure II. 1 illustrates how quantizing can be done differently for each category 
thereby enabling more accurate classification by the following table 
look-up rule; (1) quantize the measurement by the quantizing rule 
for category one (2) use the quantized measurement as an address in 
a table and test if the entry is a binary one or binary zero, (3) if it 
is a binary one assign the measurement to category one; if it is a 
binary zero, repeat the procedure for category two. 



Figure I 1 . 2 i T 1 ustrates a graph of the number of reserved decisions 

versus probability threshold a. 
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Ill Texture 

Spatial environments can be understood as being spatial distributions of 
various area-extensive objects having characteristic size and reflectance or emissive 
qualities. The spatial organization and relationships of the area-extensive objects 
appear as spatial distributions of grey tone on imagery taken of the environment. We 
call the pattern of spatial distributions of grey tone, texture. 

Figure IIl.l, taken from Lewis (1971), illustrates how texture relates to 
geomorphology. There are some plains, low hills, high hills, and mountains in 
the Panama and Columbia area taken by the Westinghouse AN/APQ 97 K-band 
radar imager system. The plains have apparent relief of 0-50 meters, the hills 
have apparent relief of 50-350 meters, and the mountains have apparent relief of 
more than 350 meters. The low hills have little dissection and are generally smooth 
convex surfaces whereas the high hills are highly dissected and have prominent ridge 
crests. 

The mountain texture is distinguishable from the hill texture on the basis of 
the extent of radar shadowing (black tonal areas). The mountains have shadowing 
over more than half the area and the hills have shadowing over less than half the 
area. The hills can be subdivided from low to high on the basis of the abruptness 
of tonal change from terrain front slope to terrain back slope. 

There have been six basic approaches to the measurement and quantifica- 
tion of image texture; autocorrelation functions (Kaizer, 1955), optical transforms, 
(Lenddris and Stanley, 1970), digital transforms, (Gramenopoulos, 1973; Hornur^g 
and Smith, 1973; Kirvida and Johnson, 1973), edgeness (Rosenfeld and Thurston, 
1971), structural elements, (Matheron, 1967; Serra, 1973), and spatial grey tone 
co-occurrence probabilities, (Haralick et al., 1973). The first three of these 
approaches are related in that they all measure spatial frequency directly or in- 
directly. Spatial frequency is related to texture because fine textures are rich in 
high spatial frequencies while coarse textures are rich in low spatial frequencies. 

An alternative to viewing texture as spatial frequency distribution is to 
view texture as amount of edge per unit area. Coarse textures have a small number 
of edges per unit area. Fine textures have a high number of edges per unit area . 

The structural element approach uses a matching procedure to detect the 
spatial regularity of shapes called structural elements in a binary image. When 
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the structural elements themselves are single resolution cells, the information 
provided by this approach is the autocorrelation function of the binary image. 

By using larger and more complex shapes, a more generalized autocorrelation can 
be computed. 

The grey tone co-occurrence approach characterizes texture by the spatial 
distribution of its grey tones. Coarse textures are those for which the distribution 
changes only slightly with distance and fine textures are those for which the dis- 
tribution changes rapidly with distance, 

lll.l Optical Processing Methods and Texture 

Edward O'Neill's (1956) article on spatial filtering introduced the engineering 
community to the fact that optical systems can perform filtering of the kind used in 
communication systems. In the case of the optical systems, however, the filtering 
is two-dimensional. The basis for the filtering capability of optical systems lies in 
the fact that the light amplitude distributions at the front and back focal planes of 
lens are Fourier Transforms of one another. The light distribution produced by the 
lens is more commonly known as the Fraunhofer diffraction pattern. Thus, optical 
methods facilitate two-dimensional frequency analysis of images. 

The paper by Cutrona et al. (1960) provides a good review of optical pro- 
cessing methods for the interested reader. More recent books by Goodman (1968), 
Preston (1972), Shulman (1970) comprehensively survey the area. 

In this section, we describe the experiments done by Lendaris and Stanley, 
Egbert et a|., and Swanlund using optical processing method? in aerial or satellite 
imagery. Lendaris and Stanley (1970) illuminated small circular sections of low 
altitude aerial photography and used the Fraunhofer diffraction pattern as features 
for identifying the sections. The circular sections represented a circular area on 
the ground of 750 feet. The major category distinction they were interested in 
making was man-made versus non man-made. They further subdivided the man-made 
category into roads, road intersections, buildings, and orchards. 

The pattern vectors they used from the diffraction pattern consisted of 40 
components. Twenty components were averages of the energy in 9° wedges of the 
diffraction pattern. They obtained over 90 percent identification accuracy. 
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Ulaby and McNaughton used an optical processing system to examine the 
texture of ERTS imagery over Kansas. They used circular areas corresponding to a 
ground diameter of about 37 km and looked at the diffraction patterns for four 
different physiographic regions in Kansas. They used a diffraction pattern sampling 
unit having 32 sector wedges and 32 annular rings to sample and measure the diffrac- 
tion patterns. (See Jensen (1973) for a description of the sampling unit and its use 
in coarse diffraction pattern analysis.) They were able to interpret the resulting 
angular orientation graphs in terms of dominant drainage patterns, roads and fields 
but interpreted the spatial frequency graphs in terms of stress patterns, rough 
terrain and field patterns. Their results indicated that the spatial frequency in- 
formation was highly correlated with physiography. 

Swanlund (1969) has done work using optical processing on aerial images 
to identify species of trees. Using imagery obtained from Itasca State Park in 
northern Minnesota, photo interpreters identified five (mixture) species of trees 
on the basis of the texture: Upland Hardwoods, Jack pine overstory/Aspen under- 
Story/Upland Hardwoods understory. Red pine overstory/Aspen understory, and 
Aspen. They achieved classification accuracy of over 90 percent. 

III. 2 Texture and Edges 

The autocorrelation function, the optical transforms, and the fast digital 
transforms (FFT and FHT) basically all reference texture to spatial frequency. 

Rosenfeid and Thurston (1971) conceive of texture not in terms of spatial frequency 
but in terms of edgeness per unit area. An edge passing through a resolution cell 
is detected by comparing the values for local properties obtained in pairs of non- 
overlapping neighborhoods boarding the resolution cell. To detect microedges, 
small neighborhoods must be used. To detect macroedges, large neighborhoods 
must be used . 

The local property which Rosenfeid and Thurston suggested was the quick 
Roberts gradient (the sum of the absolute value of the differences between diagonally 
opposite neighboring pixels). Thus, a measure of texture for any subimage is obtained 
by computing the Roberts gradient image for the subimage and from it determining the 
average value of the gradient in the subimage. Triendl (1972) uses the Laplacian 
instead of the Roberts gradient. 
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Sutfon and Hall (1972) extended Rosenfeld and Thurston's Idea by making the 
gradient a function of the distance between the pixels. Thus, for every distance d 
and subimage I defined over a neighborhood N of resolution cells, they compute 

g(d)= 2 K' + d/j)l + !K> J) 

(i/j) N 

+ |I('/]) " I('/j + d)| + jl(i,j) - I(i,j”d)|} . 

The graph of g(d) is like the graph of the minus autocorrelation function translated 
vertically. 

Sutton and Hall applied this textural measure in a pulmonary disease 
identification experiment using radiographic imagery and obtained identification 
accuracy in the 80 percentile range for discriminating between normal and abnormal 
lungs when using a 128 x 128 subimage. 

III. 3 Digital Transform Methods and Texture 

In the digital transform method of texture analysis, the digital image is 
typically divided into a set of non~overlapping small square subimages. Suppose 
the size of the subimage is n x n resolution cells, then the n^ grey tones in the 
subimage can be thought of as the n^ components of an n^~dimensionaI vector. 

In the transform technique, each of these vectors Is re~expressed In a new coor~ 
dinate system. The Fourier Transform uses the sine-cosine basis set. The Hadamard 
Transform uses the Walsh function basis set, etc. The point to the transformation 
is that the basis vectors of the new coordinate system have an interpretation that 
relates to spatial frequency (sequency) and since frequency (sequency) is a close 
relative of texture, we see that such transformation can be useful. 

Gramenopoulos (1973) used a transform technique using the sine"cosine 
basis vectors (and implemented it with the FFT algorithm) on ERT5 imagery to 
investigate the power of texture and spatial pattern to do terrain type recognition. 
He used subimages of 32 by 32 resolution cells and found that on Phoenix, Arizona 
ERTS image 1940-17324-5 spatial frequencies larger than 3.5 cycles/km and 
smaller than 5.9 cycles/km contain most of the information needed to discriminate 
between terrain types. The terrain classes were: clouds, water, desert, farms, 
mountains, urban, riverbed, and cloud shadows. He achieved an overall identifi- 
cation accuracy of 87 percent. 

Hornung and Smith (1973) have done work similar to Gramenopoulos but 
with aerial multispectral scanner imagery instead of ERTS imagery. Maurer (1974) 
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used Fourier series analysis on some color aerial film to obtain textural features to 
help determine crop types. 

Kirvida and Johnson (1973) compared the fast Fourier, Hadamard, and 
Slant Transforms for textural features on ERTS imagery over Minnesota. They used 
8x8 subimages and five categories: Hardwoods, Confiers, Open, Water, City. 
Using only spectral information, they obtained 74 percent correct identification 
accuracy. When they added textural information, they increased the identification 
accuracy to 99 percent. They found little difference between the different trans- 
form methods. 


III. 4 Spatial Grey Tone Dependence: Co-occurrence 

One aspect of texture is concerned with the spatial distribution and spatial 
dependence among the grey tones in a local area. Darling (1968) used statistics 
obtained from the nearest neighbor grey tone transition matrix to measure this 
dependence for satellite images of clouds and was able to identify cloud types on 
the basis of their texture. Read and Jayaramamurthy (1972) divided an image into 
all possible (overlapping) subimages of reasonably small and fixed size and counted 
the frequency for all the distinct grey tone patterns. This is one step more general 
than Dorling but one that requires too much memory if the grey tones can take on 
very many values. Haralick (1971) and Haralick et al . (1972, 1973) suggested an 
approach which is a compromise between the two. He measures the spatial depen- 
dence of grey tones in a co-occurrence matrix for each fixed distance and/or 
angular spatial relationship and uses statistics of the matrix as measures of image 
texture. 


t*h 

The co-occurrence matrix P = (p.j) has its (i,j) entry P.. defined as the 
number of times grey tone i and grey tone j occur in resolution cells of a subimage 
having a specified spatial relation, such as distance 1 neighbors. The textural 
features for the subimage are obtainable from the co-occurrence matrix by measures 
such as o 
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Haralick et al . (1973) list 14 different kinds of measures. 

Using statistics of the co“occurrence matrix, Haralick performed a number 
of identification experiments. On a set of aerial imagery and eight terrain classes 
(old residential, new residential, lake, swamp, marsh, urban, railroad yard, scrub 
or wooded), he obtained 82 percent correct identification with 64 x 64 subimages. 

On an ERTS Monterey Bay, California, image, he obtained 84 percent correct 
identification using 64 x 64 subimages and both spectral and textural features on 
seven terrain classes: coastal forest, woodlands, annual grasslands, urban areas, 
large irrigated fields, small irrigated fields, and water. On a set of sandstone 
photomicrographs, he obtained 89 percent correct identification on five sandstone 
classes: Dexter“L, Dexter~H, St. Peter, Upper Muddy, Gaskel , 

The wide class of images on which they found that grey ione co-occurrence 
carries much of the texture information is probably indicative of the power and 
generality of this approach. 

III. 5 A Textural Transform 

Each of the approaches described for the quantification of textural features 
had the common property that the textural features were computed for sub images of 
typical sizes such as 8 x 8, 16 x 16, 32 x 32, or 64 x 64 resolution cells. To 
determine the textural features for one pixel we would naturally center a subimage 
on the specified resolution cell and compute the textural features for the subimage. 

If we had to determine the textural features for each pixel in an image we would 
be in for a lot of computation work and would significantly increase the size of our 
data set. Thus, the usual approach has been to divide the image Into mutually 
exclusive subimages and compute the textural features on the selected subimages. 
Unfortunately, this procedure produces textural features at a coarser resolution 
than the original image. 

In this section we generalize the grey tone co-occurrence textural feature 
extractor to the textural transform mode and show how by only doubling or tripling 
the computation time required to determine the grey tone co-occurrence matrix it 
is possible to produce a resolution perserving textural transform in which each 
pixel in the transformed image has textural information about its own neighborhood 
derived from both local and global grey tone co-occurrence in the image. This 
kind of textural transform is in the class of image dependent non-| inear spatial filters. 
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Let X be the set of resolution cells of an image I (by rowcolumn 
coordinates). Let G be the set of grey tones possible to appear on image I. Then 
I: Z^ X Z^->- G. Let R be a binary relation on Z^ x Z^ pairing together all those 
resolution cells in the desired spatial re lat ion. The co-occurrence matrix 
P, P; G X G [0,1], for image I and binary relation R is defined by 

p(i j) _: ^ i((a,b),(c,d))eR|l(g,b)= i and I(c,d>=i } 

The textural transform J,J: Z x Z of image I relative to function f, is 

r c 

defined by 


J(y,x) = 


1 


E f[Pa(y,x),I(a,b))l 
(a,b)eR{y,x) 


Assuming f to be the identity function, the meaning of J(y,x) is as follows. 
The set R(y,x) is the set of all those resolution cells in Z^ x Z^ in the desired spatial 
relation to resolution cell (y,x). For any resolution cell (a,b)€R(y,x), 
P(I(y,x),I(a,b)) is the relative frequency by which the grey tone I(y,x), appearing 
at resolution cell (y,x), and the grey tone I(a,b), appearing at resolution cell 
(a,b), co-occur together in the desired spatial relation on the entire image. The 
sum 

2 P(I(ypx)/ 1(0/ b)) 

(a,b)eR(y,x) 

is just the sum of the relative frequencies of grey tone co-occurrence over all 
resolution cells in the specified relation to resolution cell (y,x). The factor 

# R(y ~ x) ' reciprocal of the number of resolution cells in the desired spatial 

relation to (y,x) is just a normalizing factor. 
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IV 


Spatial Pre-Processing 

Spatial enhancement processes can be implemented before or after the 
classification of the original images. One spatial averaging process which can 
be used before classification of the original image is rectangular convolution. 

A 2 X 2 rectangular convolution, for example, is the process that replaces the 
left upper resolution cell of each 2x2 window by the average of the grey tones 
in the 2x2 window. A 3 x 3 rectangular convolution replaces each grey tone 
with the average of the grey tones in a 3 x 3 window centered around it. The 
process of rectangular convolution can be implemented before or after texture 
transform. The window size for the rectangular convolution process can be as 
big as required. 

Figure IV illustrates how the rectangular convolution can enhance the 
textural transform processed images. Notice that the rectangular region on the 
left lower corner is not easy to distinguish on the image with no rectangular 
convolution before or after texture transform. Figure IV a, but it is distinguishable 
on Figure IVd, the image with 2x2 rectangular convolution before texture 
transform and no rectangular convolution after texture transform, as it is on 
Figures IVe to IV i. The two strips on the middle of the image are not easily 
distinguished on Figures IVa to IVf, but they are easily distinguished on Figure 
IVg, the image with 3x3 rectangular convolution before texture transform and 
no rectangular convolution after texture transform. They are also distinguishable 
on images IV h and IV i which have been processed with a 3 x 3 convolution after 
the textural transform. For distinguishing rectangular region and the two strips 
on the image. Figure IV i, the image with 3x3 rectangular convolution before 
and after texture transform seems best. 
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V 


spatial Posl— Processing 

Spatial post processing the classified image can be used to reduce image 
complexity and achieve some degree of spatial simplification and generalization. 

Two post processing techniques are region filling and shrinking. A region filling 
operation assigns an unassigned resolution cell to the category assignment of one 
of its neighboring resolution cells. 

A resolution cell can be defined to have the four resolution cells above, 
below, to the left, and to the right of it as neighbors or to have those plus the 
resolution cells diagonally neighboring it as its neighbors. The first set of 
resolution cells is called its 4“neighbors and the second set of resolution cells 
is called its 8-neighbors. The concepts of 4-neighboring and 8-neighboring is 
illustrated in Figure V.l. 

A region filling operation which assigns an unlabeled resolution cell to 
the category assignment of one of its four nearest neighbors is called a 4-fill 
operation. A region filling operation which assigns an unlabeled resolution cell 
to the category assignment of one of its eight nearest neighbors is called an 
8-fill operation. A region filling operation which iterates first filling using 4 
neighbors and then 8 neighbors then 4 then 8 etc., until all resolution cells are 
labeled, we shall for simplicity call region filling. 

Figure V.2 illustrates the advantage of region filling alternating between 
4-neighbors or 8-neighbors. A labeled resolution cell in an area of unlabeled 
resolution cells would grow as a diamond region under repetitive 4-fill operations. 

It would grow as a square region under repetitive 8-fill operations. And it would 
grow almost as a circle under repetitive 8-fill and 4-fill operations. 

Region shrinking is the opposite kind of operation from region filling. A 
region shrinking operation assigns a labeled resolution to "unassigned" if its 
neighbors have different labels from it. 

A region shrinking operation which assigns a labeled resolution cell to 
"unassigned" if k of its four nearest neighbors have labels which are different than 
its own label is called a 4~k shrink operation. A region shrinking operation which 
assigns a labeled resolution cell to "unassigned" if k of its eight nearest neighbors 
have labels which are different from its own label is called and 8-k shrink operation. 
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Figure V . 1 
Figure V . 1 


a il I usfrqtes the 4-neighborhood 
of a resolution cell and 

b illustrates the 8-neighborhood 
of a resolution cell . 
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4 and 8“fills 


Figure V.2 illustrates the effect of 4 and 8-filling 
or a single resolution cell. 
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In Figure V.3 we illustral-e the effect of the filling and shrinking operations 
on a classified image. Figure V.3 a is a classified image. The black areas represent 
uiiassigned resolution cells. (The decision rule leaves unassigned those resolution 
cells having multispectral signatures which do not provide enough information to 
make a reliable assignment.) Figure V.3b shows the classified image of Figure V.3a 
after a complete region filling. Notice that after a complete region filling, all 
resolution cells have a label. Figure V.3c shows the classified Image of Figure V.3a 
after a 4-0 shrink. Notice that it has more black area than the image in Figure V.3a 
due to the effect of its relabeling labeled resolution cells to "unassigned" . 
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' Example of a non“Convex set which Is not regular 
Figure 3 illustrates the relationship between set convexity and regularity 
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VI Spectral Anolysis: Edit 6 

Of the 6 best spectral bands on edit ^6, ,40 “ ,44, ,588 “ .643, .65 ~ 

.69, ,72 ■■ .76, .981 “ 1.045, and 2.10“ 2,36 micrometers, the feature selection 
procedure' selected band pairs .40 - ,44 and .65 “ ,60 with .40 “ ,44 and 2,10 - 
2.36 micrometers as the best 2 band pairs for the table look“up rule. Figure VI. 1 
shows the ,72 “ .76 micrometer band and Figure VI, 2 shows the ground truth 
training data overlay on this band. The alpha-beta thresholds were set at ,3 and 
.021 . This threshold selection was too low for of the 159,500 points to be 
classified, 67,323 were reserved assignments because of incompatible assignments 
between the first and second band pairs and 6,928 were reserved assignment because 
there was more than one possible assignment common to the two band pairs. Figure 
VI, 3 shows the resulting classification. The contingency table. Table VI, 1 shows 
an equally weighted mlsidentification error rate of 36% and equally weighted false 
identification error rate of 34%, The largest cause of the mlsidentification error 
was category 2.4, immature poletimber loblolly pine, being assigned to Category 
1 .3, immature sawtimber shortleaf pine, and category 2.6, mature sawtimber 
loblolly pine being assigned to category 2,5, immature sawtimber loblolly pine 
and being assigned to category 2.3, seedling and sapling loblolly pine. 

If the classified image is filled so that all resolution cells whose category 
assignment was reserved is assigned to the category of its spatially nearest resolution 
cell neighbor which is assigned, the error rate remains substantially the same, about 
a 36% mlsidentification and false identification error rate (Figure VI ,4 and Table 
VI .2). This implies that for those resolution cells whose assignment was reserved 
because of the low probability of corf* ct assignment, category assignments, 
almost as good as those originally assigned, can be made using the spatial in~ 
formation carried by the initially classified image with the reserved decisions. 

Perhaps what is even more surprising about the amount of spatial information 
the classified image has is that by performing spatial operations on it, the- classifico” 
tion accuracy can increase. For example, if the completely filled image is shrunk 
for one iteration with a simple 4-shrink operotor and then filled up again. Table VI, 3 
shows an accuracy increase: 33% mlsidentification error rate and 35% false 
identification error rate. Comparable results are also obtained by using the initially 
classified image with reserved decisions and performing a 4“fill iteration followed 
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by an 8-fiII iteration followed by a 4-shrink iteration and then completely filled 
(Figure VI. 5 and Table VI o4)o 

The best (percentage wise) 2 band pair results came from starting with the 
initially classified image with reserved decisions and doing a 4-fill, an 8-fill, a 
4-shrink, an 8-shrink, and then a complete filling upo This yields a 31% mis- 
identification error rate and 7% false identification error rate (Table VI. 5 and 
Figure VI. 6). Notice, however, that all the points in category 2.4, poletimber 
immature loblolly, have been misidentified as category 1 o3, sowtimber immature 
shortleaf pine, and all the points in category 2.6, mature sowtimber loblolly pine, 
have been misidentified as categories 1 .3, 2.3 and 2.5. Furthermore, no points 
were assigned to categories 2.4 and 2,6 o This suggests that the tree stands in 
those areas of immature loblolly and mature sawtimber loblolly pine had a sub- 
stantial number of trees spectrally similar to those in categories 1 .3, 2.3, and 
2o5. Areas predominantely in categories 2.4 and 2.6 would have some resolution 
cells initially assigned to categories 2.4 and 2.6 plus wrong assignments to 
categories 1 .3, 2.3, or 2.5. Hence, a context sensitive shrinking operation on 
the 4-fill and 8-fill image which would leave alone any resolution cell assigned 
to category 2.4 if it neighbors a resolution cell of category 1 .3 and which would 
leave alone any resolution cell assigned to category 2.6 if it neighbors a resolu- 
tion ceil of category 1.3, 2.3 or 2.5 has the possibility of permitting a higher 
probability of correct identification. 

If instead of doing only one 4-shrink then 8-shrink iterations, two such 
iterations are made before a complete filling, then the results are not quite as 
good: 34% misidentification error rate and 6% false identification error rate. 
(Table VI.6). 

The use of additional spectral bands can sometimes increase identification 
accuracy. In the case of the edit ^6 data, this did not seem to be the case. The 
three best band pairs were: 

(1) .40 - .44 and .65 - .69 micrometers 

(2) .40 - .44 and 2.10 - 2.36 micrometers 

(3) o72 - .76 and .981 - 1 .045 micrometers 

The alpha-beta thresholds were set at .6 and .042, respectively. The resulting 
number of reserved decisions due to no common category assignment was 51,794 
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and the number of reserved decisions due to more than one possible category 
assignment was 19,706 (Figure VI .7 and Table VI .7) o Higher thresholds would 
have been better. 

After a complete filling, there was a 34% misidentification and 33% false 
identification error rate (Figure VI. 8 and Table VI. 8). If the completely filled 
image had a 4-shrink operation and then another complete filling, the misidenti- 
fication error rate Improved to 31% and false identification error rate improved 
to 16% (Figure VI. 9 and Table VI«9)o If before the complete filling is done an 
iteration of a 4-fill followed by an 8”fill and a 4“shrink followed by an 8”shrink 
is done, the misidentification error rate inproves to 30% and the false identifica- 
tion error rate improves to 5%, the best 3-band pair result (Figure VI. 10 and 
Table VI.10)o As in the two band pair case, doing two iterations of the 4-shrink 
followed by the 8-shrink instead of one iteration, does not provide as much 
improvement: a 36% misidentification error rate and a 6% false identification 
error rate (Table VIoll). The best 3 band pair result confused the same categories 
as the best 2 band pair result. Category 2 .4, poletimber immature loblolly was 
assigned as category 1 .3, immature shortleaf pinCo Category 2.6, mature sow- 
timber loblolly pine was assigned to categories 2^3 and 2.5, seedling and 
sapling loblolly and sawtimber immature loblolly pine. 
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Figure VI. 1 The .72 ~ .76 micromefer band 



Figure VI. 2 


The ground truth training data overlayed on the .72 ~ .76 
micrometer band. 
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Figure VI. 3 The classification of the best two band pairs for alpha 
beta thresholds of .3 and .021 . 



Figure VI. 4 The classified image of Figure VI. 3 after a complete 
filling . 
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Table VI. 1 The contingency table of the best 2 band pairs for alpha- 
beta thresholds of .3 and .021 . 
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Table VI. 2 The contingency table of the best 2 band pairs after a complete 

filling. 
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Table VI. 3 The contingency table of the best 2 band pairs after complete 
filling, 4-shrink, and complete filling operations. 
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Table VI. 4 The contingency table of the best 2 band pairs after 4-fill, 
8-fiII, 4-shrink, and complete filling operations. 
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Figure VI. 6 The classified image of Figure VI. 3 after 4-fill, 8”fill, 
4”shrink, 8“shrink, and then complete filling operations. 
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Table VI. 5 The contingency table of the best 2 bond pairs after 4-fill, 

8-fill, 4“shrink, 8-shrink, and complete filling operations. 
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Table VI. 6 The contingency table of the best 2 band pairs after 4-fiil, 
8-fill, 4-shrink, 8-shrink, 4-shrink, 8-shrink and complete 
filling operations. 
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Figure VI, 7 The classification of the three best band pairs for alpha 
beta thresholds of .6 and .042. 



Figure VI. 8 The classified image of Figure VI. 7 after a complete 
filling . 


36 



CONTINGENCY TABLE FOR SAMH21GDT - 1 SAMH?BB02 - 1- SCALE FACTOR-IO** 0- 


COL. - ASSIGN CAT ROW_; LRuE CAX 



R DEC 

1.3 

1.4 2.3 

2.4 2 . 5 ' 

2.6 

7 . 2 ' 

TOTAL 

ERir 

err 


1 5Rq A 

93 O 2 14353 

390 L5106 

498.8 

9713 

130425 

(1 

0 

1.3 

2903 

33S7 

113 74 

66 157 

117 

- 139 

6956 

66it 

16 

l.A 

75') 

34 

1807 7 

0 43 

22 _ 

7 

2670 


6 

2.3 

2541 

139 

7 4521 

2 339 

137 

41 

7727 

665 

13 

2.4 

311 

217 

13 4 

23 . 31 

. 2 l_ 

9 . 

629_ 

_,295 

93 

2 .S 

1879 

260 

"20 eV" 

3 1535 

246 

7 

4034 

620 

29 

2.6 

754 

137 

6 83 

2 290 

152_ 

10 

1434 

528 

78 

7.2 

1683 

115 

119 47 

2 1 

5 

3653 

5625 

289 

7 

TOTAL71SOO 

201B3 

11387 19173 

488 17502 

5688 

13579 

159500 

3176 

34 

ERR 

0 

902 

278 299 

75 861 

548 

213 

3176 


»**«* 


ERR.. 0 21 13 S _ 77 _ 36 7.8 6 33 ***** 


Table VI. 7 The contingency table of the best 3 band pairs for alpha - 
beta thresholds of .6 and .042. 
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Table VI. 8 The contingency table of the best 3 band pairs after a 
complete filling. 
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Figure VI. 9 The classified image of Figure VI. 8 after a 4-shrink 
operation and then a complete filling. 



Figure VI. 10 The classified image of Figure VI. 7 after 4-fill, 8“fill, 
4-shrink, 8“shrink and complete filling operations. 
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Table VI. 9 The contingency table of the best 3 barrel pairs after 

complete filling, 4-shrink, and complete filling operations. 
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Table VI. 10 The contingency table of the best 3 band pairs after 4-fill, 
8“fill, 4“shrink, 8-shrink, and complete filling operations. 
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Table VI. 11 The confingency table of the best 3 band pairs after 4-fil|, 
8-fill, 4-shrink, 8-shrink, 4-shrink, 8-shrink and complete 
filling operations. 
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VII Spectral Analysis: Edit 9 


Using the same initial six spectral bands to select features from, the 
feature selector chose band pairs .40 - .44 and .65 - .69 with .72 - .76 and 
,981 “ 1.045 micrometers as the best 2 band pairs for the table look“up rule. 

Figure VII. 1 shows the .72 - ,76 micrometer band and Figure VII, 2 shows the 
ground truth training data overlayed on this band. The alpha-beta thresholds 
were set at .3 and .021 , 

The contingency table (Table VII. 1) for the best 2 band pairs classification 
with an alpha threshold of .3 and a beta threshold of ,021 gave a misidentification 
error rate of 22% and a false identification error rate of 32%. There were 79,670 
reserved assignments because of incompatible assignments between the first and 
second band pairs and 2,357 were reserved assignments because there was more 
than one possible assignment common to the two band pairs. The raw classified 
image is shown in Figure VU,3. The main cause of error is the confusion between 
category 1 .3, shortleaf pine, and category 2.5, loblolly pine. This error is due 
to assigning category 1 .3 when the true category is 2,5. A look at the timber 
stand map for edit ^9 shows a patch of category 2.5, which is surrounded by 
category 1 .3, in the lower right-hand corner. It is this area that gets mis- 
assigned the most. 

If the classified image is filled so that all resolution cells whose category 
assignment was reserved is assigned to the category of its spatially nearest 
resolution cell neighbor which is assigned, the error rate remains substantially 
the same, about a 25% misidentification error rate and 32% false identification 
error rate (Figure VII .4 and Table VII. 2) If we do 6 iterations of 4-fills and then 
do a 4-shrink and fill up, the resulting contingency table is Table VII, 3. The 
misidentification and false identification error rates of 21% and 26% are lower 
than before, but the misidentification error rate category 2,5 went from 43% to 
44% with category 1 .3 still the problem. 

The best 2 band pair results were obtained from doing a 4-shrink following 
the original classification and then filling (Figure VII. 5). Table VII.4 shows a 
misidentification error rate of 14% and a false identification error rate of 17%, 
but still the misidentification of category 2.5 is the main cause of error. The 
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shrinking first does eliminate a significant amount of error between category 3.1, 
laurel oak, and category 4.2, low quality sweetgum. Neither procedure has 
trouble classifying category 2.5 on the left-side of the timber stand. Only on 
the right side where category 2.5 resembles category 1 .3 spectrally is there 
confusion. This confusion could be ultimately due to sun angle. 

The three best band pairs were: 

(1) .40 - .44 and .65 - .69 micrometers 

(2) .72 - .76 and .981 - 1 .045 micrometers 

(3) .40 - .44 and 2.10 - 2.36 micrometers 

Figure VII. 6 shows a plot of the alpha threshold versus the number of reserved 
decisions. For the three best band pairs, the alpha and beta thresholds that 
minimized the number of reserved decisions was .6 and .042, respectively. 

The raw classified image is shown in Figure VII. 7. The contingency table 
indicates a misidentification error of 24% and a false identification error of 
30% (Table VII. 5). 

After a complete filling, there was a 25% misidentification and 32% 
false identification error rate (Figure VII. 8 and Table VII. 6). If instead, our 
post processing consisted of a 4-fill, 8-fill, 4-shrlnk, 8-shrink and then a 
complete filling the misidentification error rate was 9% and the false- identi- 
fication error rate was 9% (Table VII. 7 and Figure VII.9). 
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Figure VII. 1 The .72 - .76 micrometer band. 



Figure VII. 2 The ground iruth training data overlayed on the .72 - .76 
micrometer bond. 
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Table VII. 1 The contingency table of the best 2 bond pairs for alpha 
beta thresholds of .3 and .021 . 
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Table VII. 2 The contiiigency table of the best 2 band pairs after a 
complete filling. 
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Figure VII. 4 The classified image of Figure VII. 3 affer a complefe 
filling . 
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Table VII. 3 The contingency table of the best 2 band pairs after complete 
tilling, 4“shrink, and complete filling operations. 
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Table VII. 4 


The contingency table of the best 2 band pairs after 4-shrink, 
and complete filling operations. 
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Figure VII. 


The classified image of Figure VII. 3 after 4-shrink and 
complete filling operations. 
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Figure VII.6 A plot of the alpha thresholds versus 
number of reserved decisions 
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Table VII. 5 
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The contingency table of the best 3 band pairs for alpha - 
beta thresholds of .6 and .042. 
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Table V1I.6 The contingency table of the best 3 bond pairs after a complete 
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Figure VII. 7 The classification of the three best band pairs for alpha 
beta thresholds of .6 and .042. 



Figure VII. 8 The classified image of Figure VII. 7 after a complete 
filling . 
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Table VII. 7 The contingency table of the best 3 bond pairs after 4-fill, 

8“fill, 4“shrink, 8-shrink, and complete filling operations. 
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Figure VII. 9 The classified image of Figure VII. 7 after 4-fill, 8-fill, 
4-shrink, 8-shrink, and complete filling operotions. 
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VIII Spectral Analysis; Edit 14 

The same six spectral bands were chosen from edit ^14 as were taken from 
edit ^6 and edit ^ 9 . Figure VIII. 1 shows the .72 ~ .76 micrometer band for edit 
14 and Figure VIII. 2 shows the selected ground truth training data. The selection 
procedure chose .40“ .44 and 2. 10 ~ 2.36 with .588- .643 and 2. 10 - 2.36 
micrometers as the best 2 band pairs for the table look-up rule. The alpha and 
beta thresholds were set at .3 and .021 respectively. The thresholds were too low 
and resulted in 56,320 reserved decisions in the contingency table for classifica- 
tion (Table VIII. 1). The resulting misidentification error rate was 28% and false 
identification error rate was 29%. The result on the best 2 band pairs with 4-fill, 
8-fill, 4-shrink, 8-shrink, and complete filling operations (Table VIII. 2), was a 
misidentification error rate of 15% and a false identification error rate of 17%. 

The feature selection procedure chose band pairs .40 - .44 and .65 - .69 
micrometers, along with the best 2 band pairs for the best 3 band pairs. Using 
alpha and beta thresholds of .6 and .042, respectively, the number of reserved 
decisions was 43,236, with 25,794 points reserved because no assignment was 
possible and 17/442 reserved due to possible multiple assignments. 

The largest cause of error for best 3 band pairs (Table VIII. 3) was the 
confusion between categories 2.3 and 2.5, different ages of loblolly pine, and 
the confusion of each of these with category 4.1, low quality sweetgum. The 
misidentification and false identification error rates (46% and 48%) for category 
4.1 are high but the number of points whose true category is 4.1 is small. Figure 
VIII. 3 shows the resulting classification. There was such a small area of sweetgum, 
category 4.1, on the timber stand map that the ground truth may not be adequate 
to allow good spectral estimation. 

The first post processing procedure we used was a complete filling (Table 
VIII. 4 and Figure VIII. 4). The errors were increased by the procedure, so one 
4“shrink operation was performed on the image and this reduced the misidentifica- 
tion error to 9% and false identification error to 4% (Table VIII. 5 and Figure VIII. 

5), but the low error rates were helped by the fact that there were 84,828 
reserved decisions. Table VlII.5 does show that the confusion with category 4.1, was 
almost eliminated, though the misidentification error rate caused by assigning 
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2.3 to 2o5, 21% was still higho Completely filling the image resulted in a mis“ 
identification error rate of 17% and false identification error rate of 13% (Table 
VIII. 6 and Figure VIII. 6). 

If on the raw classified image we do one 4”fill (Table VIIlo7 and Figure 
VIII o7) and then one 8-filt, the resulting contingency table (Table VIII o8 and 
Figure VIII .8) is almost identical to Table VIII a3o The error rates on each are 
exactly the same. Then doing a 4-shrink (Table VIII. 9 and Figure VIII o9) we find 
a contingency table almost identical to Table VIII .4o But if instead of filling 
we do an 8“shrink, we almost totally eliminate error (Table VIII olO and Figure 
VIII, 10). Only 2 points are incorrectly identified. Now if we completely fill 
the image we get our best results (Table VIII, 11 and Figure Vlll.ll): 13% 
misidentification and 9% false identification error rates „ Visual comparisons 
show the closeness of the two operations « Following the fills with a 4-shrink 
produces Figure Vlll.S. Figure VIII „6 is the final classified image after complete 
filling, a 4“shrink and then a complete filling, while Figure VIII oil is tl?e final 
result of a 4-fill, 8-fill, 4-shrink, 8-shrink, and complete filling. From the 
figures, we can see that the extra shrink allowed the categories to be more dense. 
The contingency table of the image should show better results since the categories 
on the timber stand map tend to be dense, which is the case. 

The results of the shrinking operations indicate that the errors that did 
occur were sparse enough to be wiped out with the shrinking. The reason that 
a shrink operation is not performed first on the image is that It tends to eliminate 
small area categories, even though correctly assigned, on the image. 


54 




Figure VIII. 1 The .72 " .76 micrometer band. 



Figure VIII. 2 The ground truth training data overlayed on the .72 ~ .76 
micrometer band. 
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Table VIII. 1 The contingency table of the best 2 band pairs for alpha ~ 
beta thresholds of .3 and .021. 
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Table VIII. 2 The contingency table of the best 2 band pairs after 4-fill, 
8-fiU, 4-shrink, 8-shrink, and complete filling operations. 
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Table VIII. 3 The confingency table of the best 3 band pairs for alpha - 

beta thresholds of .6 and .042. 
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Table VIII. 4 The contingency table of the best 3 band pairs after a complete 

filling . ^ 





Figure VIII. 3 The classification of the three best band pairs for alpha- 
beta thresholds of .6 and .042. 
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Table VIII. 5 The contingency fable of the best 3 bond pairs after complete 
filling and 4"shrink operations. 


contingency table for SAMH4 GOT - 1 SMK4F2B06 - 1 SCALE FACTOR 10** 0 
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Table VIII. 6 The contingency table of the best 3 band pairs after complete 
filling, 4-shrink, and complete filling operations. 
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Figure VIII. 5 The classified image of Figure VIII. 4 after a 4-shrink 
operation. 
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Table VIII. 7 The contingency table of the best 3 band pairs after a 4-fill 
operation . 
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Table VIII. 8 The contingency table of the best 3 bond pairs after 4-fill, 
and 8-fill operations. 

ORIGINAL PAGH'lg 
OF POOR QUALITY 
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Figure VIII. 8 The classified image of Figure VIII, 3 after 4-fill 
and 8-fill operations. 
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Table VIII. 9 The^ contingency table of the best 3 band pairs after 4-fil| 
8-fill and 4-shrink operations. ' 
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Table VIII. 10 The^ contingency table of the best 3 band pairs after 4-fill, 

8-fill, 4-shrink, and 8-shrink operations. 
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Figure VIII. 9 The classified image of Figure VIII. 3 after 4-fill, 
8“fill, and 4-shrink operations. 



Figure VIII. 10 The classified image of Figure VIII. 3 after 4-fill, 
8-fill, 4-shrink and 8“shrink operations. 
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Table VIII. 11 The conHngency table of the best 3 band pairs after 4-fill, 
8”fill, 4-shrink, 8-shrink and complete filling operations. 
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IX Speci^|•al Analysis: Edit 3 

As with the other edits, the same 6 spectral bands were chosen, .40 ~ .44, 
.588 " .643, .65 ~ .69, .72 - .76, .981 ~ 1 .045, and 2.10 - 2.36 micrometers. 
Figure IX. 1 shows the .72 - .76 micrometer band of edit 3 and Figure IX. 2 shows 
the selected ground truth training data. 

The feature extractor chose bands .40 ~ .44 and .588 - .643 with .588 - 
.643 and .65 - .69 micrometers as the best 2 band pairs. To minimize the total 
number of reserved decisions and to try and equalize the number of reserved 
decisions due to more than one assignment and no assignment, classification for 
the two best band pairs was done using a variety of alpha and beta thresholds. 
Figure IX. 3 is a graph of the thresholds versus the number of reserved decisions. 

Table IX. 1 is the contingency table for best 2 band pairs with .3 and .021 
alpha and beta thresholds, respectively. The resulting error rates of 36% mis~ 
identification and 38% false identification are better than the corresponding 
error rates of 37% and 41% for the classification with alpha, beta thresholds 
of .4, .028 (Table IX. 2) and the corresponding error rates of 37% and 40% for 
the classification with alpha, beta thresholds of .5, .035 (Table IX. 3). But 
the total number of reserved decisions for the .3 and .021 thresholds is 
47,749. This is the highest number of reserved decisions and the lower error 
rates could be caused by lack of assignments. In this case, the fill operations 
would tend to propagate the error. Therefore, we chose .5 and .035 thresholds 
to work with. The raw classified image was post processed with 4-fili, 8-fill, 
4“shrink, 8~shrink and complete filling operations. The resulting contingency 
table (Table IX. 4) indicates an 18% misidentification error and 27% false 
identification error. The major confusion was poletimber immature shortleaf 
pine being classified as sawtimber immature shortleaf pine or poletimber immature 
loblolly pine. 

The three best band pairs consisted of the two best band pairs plus band 
pair .40 - .44 and .65 - .69 micrometers. To minimize the total number of 
reserved decisions and to try to equalize the two causes for reserved decisions, 
classification was done for the three best band pairs using a variety of alpha beta 
thresholds. The resulting graph (Figure IX, 4) indicates good alpha beta thresholds 
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are .5 and .035. Contingency table (Table IX. 5) shows a 34% misidentification 
rate and 38% false identification rate with 48,475 reserved decisions. Figure 
IX. 5 shows the resulting classification. Category 1.2 was the largest cause of 
error. It was confused with category 1 .3, sawtimber Immature shortleaf pine 
and categories 2.4 and 2.6, two kinds of loblolly pine. 

A 4-fill and an 8-fill operation reduces the misidentification error rote 
but propagates the false identification error rate (Table IX. 6 and Figure IX. 6). 
Doing a 4-shrink reduces the error rates to 18% and 23% for misidentification 
and false identification. This is as expected since fewer assignments are made 
to spatially uncertain categories but the misidentification error rate for category 
2.1 was not reduced (Table IX. 7 and Figure IX. 7). The final 8“shrink and 
then fill all the way up results in a misidentification error rate of 14% and a 
false identification error rate of 25% (Table IX. 8 and Figure IX. 8). Most of 
the error is due to category 1 .2 being confused with categories 1 .3, 2.4, and 
2.6. Thus, category 1 .2 has a misidentification error rate of 60% compared to 
6% for the next most highly confused category. Most of the confusion is 
between subclasses in the same class rather than between classes. Contingency 
table IX. 9 shows the resulting classification when categories 1.2 and 1.3 are 
combined and categories 2.4 and 2.6 are combined. The misidentification 
error rate is 10% and the false identification error rate is 14%. 
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Figure IX. 1 The .72 - .76 micromefer band. 



Figure IX. 2 The ground fruth training data overlayed on the .72 - .76 

micrometer bond. 
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Number of Reserved Decisions 



Figure IX . 3 Number of reserved decisions os a function of^ 

probability threshold alpha for best 2 band pairs, 
spectral only for edit "’o 
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Table IX. 1 The contingency table of the best 2 band pairs for alpha “ 
beta thresholds of .3 and .021 . 
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Table IX. 2 The contingency table of the best 2 band pairs for alpha - 
beta thresholds of .4 and ,028. 
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Table IX. 3 


The contingency table of the best 2 band pairs for alpha - 
beta thresholds of .5 and *035. 


CONTINGENCY TABLE FOR SAMHl GOT - 1 SMH1F3B01 - 1 SCALE FACTOR 10-»* C 


COL. - ASSIGN CAT ROW - TRUE CAT 


R 

DEC 

1. 2 

1. 3 

2. 4 

2. 6 

7.1 

TOTAL 

«ERR 

X ERR X 

SD 

UNKWN 

0 1 

1586 

8018 

25376 

3425 

21950 

70355 

0 

0 

0 

1.2 

0 

604-1 

2853 

1727 

612 

835 

12068 

6027 

50 

0 

1. 3 

0 

44 

1760 

59 

0 

0 

1863 

103 

6 

0 

2. A 

0 

0 

0 

9476 

0 

496 

9972 

496 

5 

0 

2. 6 

0 

968 

0 

269 

4632 

536 

6405 

1773 

28 

0 


7.1 

0 

0 

0 

124 

0 

4045 

4169 

124 

3 

0 

TOTAL 

0 

18639 

12631 

37031 

8669 

27862 

104832 

8523 

18 

0 

«ERR 

0 

1012 

2853 

2179 

612 

1867 

8523 

***♦♦ 

***** 

***** 

X ERR 

0 

14 

62 

19 

12 

32 

27 

-****« 

««**« 



Table IX. 4 


The contingency table of the best 2 band pairs for alpha - 
beta thresholds of .5 and .035 after a complete filling . 
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spectral only for edit ^3 
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Table IX. 5 


The contingency table of the best 3 band pairs for alpha - 
beta thresholds of .5 and .035. 
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Table IX .6 The contingency table of the best 3 band pairs after 4-fi|| 
and 8“fil| operations. 
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Figure IX. 5 The classification of the three best band pairs for alpha 
beta thresholds of .5 and .035. 



Figure IX. 6 The classified image of Figure IX. 5 after 4-fill and 8-fill 
operations. 


75 
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Table IX. 7 The contingency table of the best 3 band pairs after 4~fill, 

S-fill and 4-shrink operations. 


CONTINGENCY TABLE FOR SAMHl ODT 
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Table IX. 8 


The contingency table of the best 3 band pairs after 4-fiM, 
8-fill, 4“shrink, 8-shrink, and complete filling operations. 
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Figure IX. 8 The classified image of Figure IX. 5 after 4-fill, 8-fill, 

4-shrink, 8-shrink and complete filling operations. 
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Col. = Assign Cat. 


Row = True Cat. 
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Spectral-Text'urql Analysis; Edit 6 

We began the spectral-textural analysis of the edit ^6 data by using five 
spectral bands and two texture bands and letting the feature selection procedure 
pick the best two and best three band pairs for the table look~up decision rule. 

The five spectral bands were: 

.40 “ .44 micrometers 
.65 - .69 micrometers 
.72 - .76 micrometers 
.981 - 1.045 micrometers 
2.10 - 2.36 micrometers 

The textural transform was done on a 3x3 convolution of the .82 - .88 micrometer 
band. A second textural information band was created by doing a 3x3 convolution 
of the initial textural transform image. 

The feature selection procedure selected the two best band pairs consisting 
of: 

(1) .40 - .44 micrometer band with the 3x3 convolution before 
and after the textural transform of the .82 - .88 micrometer 
band 

(2) .65 - .69 and .981 - 1.045 micrometer bands. 

The alpha-beta thresholds were set at .3 and .021, respectively. This threshold 
selection was too low for of the 159,500 points to be classified, 74,326 were 
reserved assignments because of incompatible category assignments between the 
first and second band pairs and 1,904 were reserved assignment because there 
was more than one possible assignment common to the two band pairs. The 
resulting contingency table, (Table X.1 and Figure X.l) shows a misidentification 
error rate of 36% and a false identification error rate of 37%. After filling the 
classified image to remove all reserved assignments, the misidentification error 
rate was 38% and false identification error rate was 39%, Table X.2 and Figure 
X.2. This is worse than the best two band pair spectral results indicating that 
either the alpha-beta thresholds used created such a high number of reserved 
decisions that the classification accuracy was lowered or that a feature selection 
procedure which minimizes a lower bound on the error rate does not necessarily 
produce the features of the best classification . 
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Spatial processing can improve the identification accuracy of the 
initially classified image. For example, if the completely filled image is shrunk 
for one iteration with a 4“shrink operator and then filled again, the misidentifi- 
cation and false identification error rates improve to 33%, Table X.3 and Figure 
X.3. The biggest cause of errors was category 2.4 being assigned to category 1 .3 
and category 2.6 being assigned to categories 1 .3, 2.3 and 2.5. A still 
greater increase in identification accuracy results if the initially classified image 
with reserved decisions is operated on with a 4-fill, then 8“fill, then 4-shrink, 
then 8-shrink operations and then filled up completely (Figure X.4). The resulting 
contingency table. Table X.4, shows a 32% misidentification error rate and 7% 
false identification error rate. This is about the same as the best two-band 
spectral results. 

Doing two iterations of a 4-shrink followed by an 8"shrink (Figure X.5) 
instead of just one iteration as described for the previous classification produces 
not as good results. Table X.5 shows a 34% misidentification error rate and 7% 
false identification error rate. 

Repeating the 2 band experiment with an alpha threshold of .5 and a 
beta threshold of .035 reduces the number of reserved decisions to 42,226 
with 25,173 reserved decisions due to no assignment and 17,053 reserved 
decisions due to multiple assignments. The resulting classification (Table X.6 
and Figure X.6) gives a misidentification error rate of 37% and a false identi- 
fication error rate of 38%. 

A complete filling of the image (Table X,7 and Figure X.7) gives a 
misidentification error rate of 38% and 39%. The main cause of error is 
assigning category 1.3 when the true category Is 2.4 and assigning 2.5 when 
the true category is 2.6. If we do a 4-shrink on the filled image and then 
completely fill it again (Table X.8 and Figure X.8) we get a misidentification 
error rate of 32% and a false identification error rate of 36%, but now categories 
2.4 and 2.6 are completely misidentified. If instead we do a 4-fill, 8-fill, 
4-shrink, 8-shrink and then completely fill up the raw classification (Table 
X.9 and Figure X.9) we get a misidentification error rate of 30% and a false 
identification error rate of only 5%. This improvement over the (,3 and 
.021) result is due to better thresholding. So, even though the raw classifica- 
tion using an alpha threshold of .3 was a few percentage points better than the 
raw classification using an alpha threshold of .5, the large number of reserved 
decisions hindered classification accuracy with the fill and shrink operations. 
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We also did a 4-filI, 8~fili, 4~shrink and complete filling (Table X.IO 
and Figure X.10) on the raw classification using alpha threshold of .5 to see if 
we were doing too much shrinking. The resulting misidentification error rate 
of 32% and false identification error rate of 36% indicates that we were not. 

The best 3 band pairs results did significantly increase the accuracy 
over the two best spectral band pair accuracy and the two best spectral “textural 
band pair results. The band pairs selected by the feature selection procedure were: 

(1) .40 ” .44 micrometer band with the 3x3 convolution 
before end after the textural transform of the .82 - .88 
micrometer band 

(2) .65 “ .69 and 2.10 ~ 2.36 micrometer bands 

(3) .72 “ .76 and .981 ~ 1 .045 micrometer bands. 

The alpha-beta thresholds were set at .7 and .049, respectively. This resulted 
in 25,590 reserved decisions due to no common category assignment and 43,889 
reserved decisions because of more than one possible category assignment. The 
thresholds were set just a little too high. 

The contingency table of the Initially classified image with reserved 
decisions is shown in Table X.ll. It indicates a 35% misidentification error 
rate and 37% false identification error rate. Completely filling the initiolly 
classified image with reserved decisions yields a misidentification error rate 
of 38% and false identification error rate of 37%. This identification accuracy 
(Table X. 12) is just below the best 3 band pair spectral results. 

If the completely filled image is operated on with one iteration of a 
4-shrink operation and then completely filled, the misidentification error rate 
improves to 29% and false identification error rate Improves to 30% (Table X.13 
and Figure X.ll). The results indicate that almost all resolution cells originally 
assigned to category 2.4 were neighboring resolution cells of a different category. 
Hence, the 4-shrink operation eliminated most of the assignments to category 2.4. 

The basically scattered assignments to category 2,4 was manifest in the 
next experiment in which we did a 4-fill, then an 8-fill, then a 4-shrink, then 
an 8-shrink and a complete filling of the initially classified image with reserved 
decisions. The contingency fable (Table X. 14 and Figure X.12) shows a 23% 
misidentification error rate and a 6% false identification error rate. These results 
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are definitely better than the corresponding three best spectral band pair results. 

The main reason for the identification accuracy increase is that most of category 2.6 
was assigned to category 2.6; only some of category 2.6 was assigned to category 
2.5 and hardly any at all to category 1.3. All of category 2.4, however, was 
misidentified as category 1.3, 

Following the pattern of the previous results, if a double 4~shrink and 
then 8“shrink operation is applied instead of a single 4-shrink and then 8-shrink, 
the classification results are not quite as good: a 39% misidentification error 
rate and 12% false identification error rate. As shown in Table X.15, category 
2.4 is misidentified as category 1 .3 and category 2.6 is misidentified as category 
2,3 and category 2.5. 
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Table X.2 


The contingency table of the best 2 band pairs for alpha 
beta thresholds of .3 and .021 after a complete filling. 
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Figure X.l The classification of the best 2 band pairs for alpha 
beta thresholds of .3 and ,021 . 


Figure X.2 The classified imqge of Figure X.l after a complete filling 
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Table X.3 The contingency table of the best 2 band pairs for alpha - 
beta thresholds of .3 and .021 after complete filling, 
4-shrink, and complete filling operations. 
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Table X.4 The contingency table of the best 2 band pairs for alpha - 

beta thresholds of .3 and .021 after 4-fill, 8-fill, 4-shrink, 
8-shrink, and complete filling operations. 
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The classified image of Figure X.l after complete filling 
4“shrink, and complete filling operations. 
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The contingency table of the best 2 band pairs for alpha - 
beta thresholds of .Sand .021 after 4-fill, 8-fill, 4-shrInk, 
8“shrinl<, 4-shrink, 8-shrink, and complete filling operations 
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Table X.6 The contingency table of the best 2 band pairs for alpha - 
beta thresholds of .5 and .035. 
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Figure X.5 The classified image of Figure X.l after 4-fill, 8-fill, 
4-shrink, 8-shrink, 4-shrink, 8-shrink and complete 
filling operations. 



Figure X.6 


The classifcation of the best 2 band pairs for alpha - 
beta thresholds of .5 and .035. 
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Table X.7 The contingency table of the best 2 band pairs for alpha - 
beta thresholds of .5 and .035 after a complete filling. 
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Table X.8 


The contingency table of the best 2 band pairs for alpha - 
beta thresholds of .5 and .035 after complete filling, 
4-shrlnk, and complete filling operations. 
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Figure X.8 The classified image of Figure X.6 after complete filling, 
4“shrink and complete filling operations. 
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Table X.9 


The contingency table of the best 2 band c 

lairs 

for a 


beta thresholds of .5 and .035 after 4-f1ir, 8-fill, 4-shrink, 
8“shrink, and complete filling operations. 
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Figure X.9 The classified image of Figure X.6 after 4-fill, 8-fill, 
4-shrink, 8-shrink, and complete filling operations. 



Figure X.IO The classified image of Figure X.6 after 4-fill, 8-fill, 
4-shrink and complete filling operations. 
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Table X. 11 The contingency table of the best 3 band pairs for alpha - 
beta thresholds of .7 and .049. 
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Table X.12 The contingency table of the best 3 band pairs after a 
complete filling * 
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Table X.13 


The^ conHngency table of the best 3 band pairs after complete 
filling, 4”shrink, and complete filling operations. 
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Table X . 

14 

The contingency table of the best 3 bond pairs 

1 after 4-fill, 




8-flll, 4-shrink, 

8-shrink and complete filling operations. 



Figure X, 11 


The classification of the best 3 band pairs for alpha - beta 
thresholds of .7 and .049 after complete filling, 4-shrink, 
and complete filling operations. 



Figure X.^2 The classification of the best 3 band pairs for alpha - beta 
thresholds of .7 and .049 after 4-fill, 8-fill, 4-shrink, 
8-shrink, and complete filling operations. 
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Table X. 15 The confingency table of the best 3 band pairs after 4"fiil, 
8-fill, 4-shrink, 8-shrink, 4-shrink, 8-shrink, and complete 
filling operations. 
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XI 


Specfral-Textural Analysis; Edit 9 

Wii-h this edit we experimented to find the best texture transforms. The 
.82 ~ .88 micrometer band was chosen as the band having the most spatial in- 
formation (Figure XI. 1). Figure XI. 2 is a 2x2 rectangular convolution of the 
.82 - .88 micrometer band and Figure XI. 3 is a 3x3 rectangular convolution of 
the band. Each of these bands were used as inputs into the texture transform. 

The resulting textural transform images are shown in Figures XI. 4, XI. 5 and 
XI. 6. Each of these were convoluted with a 2x2 window size (shown in Figures 
XI. 7, XI. 8, XI. 9). Finally the textured transforms were convoluted with a 
3x3 convolution window giving us 3 more texture images (Figures XI. 10, XI. 11 
XI. 12). Using our own visual discretion we chose the textural transform with a 
3x3 rectangular convolution after and the 3x3 rectangular convolution before 
transforming with a 3x3 rectangular convolution after transforming as the two 
texture bands with the most information (these are shown in Figures XI. 10 
and XI. 12). 

We combined these 2 texture bands with the spectral bands and the 
feature selector chose band pairs .40 - .44 micrometers and the 3x3 rectangular 
convolution before and after the textured transform with .65 - .69 and 2.10 “ 

2.36 micrometers as the 2 best band pairs for classification. Band pair .72 - .76 
and .981 " 1 ,045 micrometers was selected with the other two for the best 3 
band pairs. Figure XI. 13 and XI. 14 show the graphs of the threshold alpha 
against the number of reserved decisions. For best 3 band pairs the best alpha 
threshold was .7 with a beta threshold of ,049. 

To check the choice of thresholds we checked several results using different 
thresholds. The best 3 band pairs classification with alpha, beta thresholds of .3 
and .021 gave us a misidentification error rate of 20% and a false identification 
error rate of 20% (Table XI, 1 and Figure XI. 15). The error rate was low but the 
total number of reserved decisions 104,531 is high. Only 89 of these points were 
reserved due to more than one assignment, while 104,443 points were reserved 
because of no assignment . The largest cause of error was due to misidentification 
of category 2.6 as category 2.5, both subclasses of loblolly pine. 
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Post processing with a 4-shrink and then a complete filling we obtained 
misidentification and false identification error rates of 36% and 20%. Both 
category 2.6 and category 3.1, laurel oak, had misidentification error rates 
of 100% (Table XI. 2 and Figure XI. 16). Though the shrink operation usually 
reduces error, if a sparse category is assigned correctly, the shrink operation 
here tended to wipe out the category. Table XI. 2 shows us that this happened to 
category 2.6 and category 3.1 . If instead of a shrink we first did a 4-fill, then 
a 4-shrink and then a complete filling, the resulting contingency table is Table 
XI. 3 (Figure XI. 17). The misidentification error rate was 18% and the false 
identification error rate was 16%, but the misidentification error rate for 
category 2.6 was still high at 41%, The main cause of error is the confusion of 
2.6 and 2.5. The only way left to eliminate the confusion is to change thres- 
holds. 

Values of .6 and .042 for the alpha, beta thresholds resulted in a mis- 
identification error rate of 25% and a false identification error rate of 28% 

(Table XI. 4). The misidentification error rate for categories 2,5 and 2.6 were 
31% and 34%, respectively. If .7 and .049 are chosen for the alpha and beta 
thresholds we get error rates of 25% and 31%, but the misidentification error 
rate for category 2.6 is only 24% and the misidentification error rate of category 
2.5 is 31% (Table XI. 5). The number of reserved decisions is 71,919 with 
43,045 points being reserved because of more than one assignment and 28,874 
points reserved because of no assignment. With thresholds for alpha and beta 
of .8 and ,063, the misidentification and false identification error rates were 
28% and 32%, respectively (Table XI. 6). Though the misidentification error 
rate for category 2.6 has been reduced to 19% and for category 2.5 it was 
reduced to 21%, the misidentification and false Identification error rates for 
category 3.1 have grown to 62% and 62%, and for category 4.2 the rates have 
gone up to 52% and 45%. In addition the number of reserved decisions has 
risen to 121,716 indicating that the thresholds hove gotten too high. 

Since the error rates for Table XI. 4 and Table XI. 5 were almost the 
same, the results from the classification with thresholds of .7 and .049 should 
be better for post processing. The main cause of error had been with categories 
2.5 and 2.5 and this classification showed lower error rates for these categories. 
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If we fill up the image with alternating 4“fill and 8‘fills we get a misidentification 
error rate of 27% and a false identification error rate of 33% (Table XI. 7). This is 
no improvement on the raw classification so the shrink operation is needed to 
eliminate incorrect assignments. Post processing with a 4-fill and an 8-fill so the 
shrink operations do not wipe out sparsely populated categories, then doing a 
4-shrink and 8-shrink and finally a complete filling, we obtain a misidentification 
error rate of 8% and a false identification error rate of 11% (Table XI. 8 and Figure 
XI. 18). The misidentification error rate for category 2.6 was reduced to 0 and 
the confusion between category 3.1 and 4.2 was small. As was the case with the 
spectral analysis ths misidentification of category 2.5 with 1.3 is the main cause 
of error. Though the texture analysis gives better overall results, it cannot over” 
come the inability of the decision rise to separate categories 2.5 and 1.3 in the 
lower right hand corner of the timber stand map. 

The results of the best 2 band pairs classification were not as good. The 
contingency table resulting from alpha, beta thresholds of .3 and .021 resulted 
in a misidentification error rate of 25% and a false identification error rate of 
31% (Table XI. 9 and Figure XI. 19). If we do a 4-fill, 4-shrink and fill up we 
get error rates of 23% and 29% (Table XI. 10 and Figure XI. 20). If we shrink 
first and then fill up, the results showed improvement with a misidentification 
error rate of 15% and a false identification error rate of 20% (Table XI. 11 and 
Figure XI. 21). 
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Figure XI. 1 The .82 


88 micrometer band used for the texture transform 


Figure XI. 2 Shows Figure XI. 1 after a 2x2 rectangular convolution 




Figure XI. 3 


Shows Figure XI. 1 after a 3x3 rectangular convolution. 



Figure XI. 4 The texture transform of Figure XI. 1. 
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Figure XI. 7 


Shows Figure XI. 4 after a 2x2 rectangular convolution. 



Figure XI. 8 Shows Figure XI. 5 after a 2x2 rectangular convolution. 
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Figure XI. 9 


Shows Figure XI. 6 after a 2x2 rectangular convolution. 



Figure XI. 10 


Shows Figure XI. 4 after a 3x3 rectangular convolution. 
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Figure XI. 11 


Shows Figure XI. 5 after a 3x3 rectangular convolution. 



Figure XI. 12 Shows Figure XI. 6 after a 3x3 rectangular convolution. 
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Figure XI. 13 Number of reserved decisions as a function of 

probability threshold alpha for best 2 band pairs 
with texture for edit ^9 
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Figure XI. 14- Number of reserved decisions as a function of 

probability threshold alpha for best 3 band pairs 
spectral only for edit ^9 
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Table XI, 1 The contingency table of the best 3 band pairs for alpha 
beta thresholds of .3 and .021 . 
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Table XI.2 The contingency table of the best 3 band pairs for alpha - 
beta thresholds of .3 and .021 after 4“shrink and complete 
filling operations. 
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Figure XI, 16 The classified image of Figure XI. 15 after 4-shrink and 
complete filling operations. 
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Table XI. 3 The conf-ingency table of the best 3 band pairs for alpha - 
beta thresholds of .3 and .021 after 4-fill, 4-shrink, and 
complete filling operations. 
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— Table XI. 4 The contingency table of the best 3 band pairs for alpha - — 
beta thresholds of .6 and .042. 


110 


Page is 
POOR QUALITY 








Figure XI. 17 


The classified image of Figure XI. 15 after 4-fill, 4-shrink 
and complete filling operations. 



Figure XI. 18 The classification of the best 3 band pairs for alpha - beta 
thresholds of .7 and .049 after 4-fill, 8-fill, 4-shrink, 

8' shrink, and complete filling operations. 
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Table XI.5 The contingency table of the best 3 band pairs for alpha 
beta thresholds of .7 and .049. 
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Table XI. 6 The contingency table of the best 3 band pairs for alpha 

beta thresholds of .Sand .063. 
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Table XI. 7 The contingency table of the best 3 band pairs for alpha “ 

beta thresholds of .7 and .049 after a complete filling. 
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Table XI. 8 The contingency table of the best 3 band pairs for alpha - 

beta thresholds of .7 and .049 after 4-fill, 8-fill, 4-shrink, - 
8-shrInk, and complete filling operations. 
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“ Table XI. 9 


The contingency table of the best 2 band pairs for alpha - 
beta thresholds of ,3 and .021 . 
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Table XI. 10 The contingency table of the best 2 band pairs after 4“fill, ~ 
4-shrink, and complete filling operations. 
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Figure XI. 19 The classification of the best 2 band pairs for alpha - 
beta thresholds of .3 and .021 . 



Figure XI. 20 The classified im<^e of Figure XI. 19 after 4-fill, 4-shrink, 
and complete filling operations. 
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- Table XI. 11 


The contingency table of the best 2 band pairs after 4-shrink 
and complete filling operations. 
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Figure XI. 21 


The classified image of Figure XI. 19 after 4-shrink and 
complete filling operations 


117 



XII Spectral-Texfural Analysis; Edit 14 

The spectral texture analysis of edit ^14 began just like that of the other 
edits except the feature selection did not choose a texture band as one of the best 
2 or best 3 band pairs. We, nevertheless, did an experiment with 2 band pairs. 

We chose bands .40 - .44 and .72 “ .76 with bands .72 - .76 and the 
texture transform image. The texture image was the result of a 3x3 convolution 
of the .82 - .88 micrometer band as input into the texture transform and a 3x3 
convolution after the texture transform. The alpha and beta thresholds were .3 
and .021, respectively. Figure XII. 1 shows the .82 ~ .88 micrometer band used 
for the texture transform. The texture transformed image that was used for pro- 
cessing is shown in Figure XII. 2. Figure XII. 3 shows the texture transform 
result with no convolution before transforming and with a 3x3 convolution after. 
The feature selector did not choose this band and visually we can see that it 
has much less spatial information than the texture transform that was chosen. 

The contingency table that resulted from the table look-up rule (Table 
XII. 1) shows a 43% misidentification error rate and a 44% false identification 
error rate. This is not nearly as good as the spectral results. There were a large 
number of reserved decisions, 72,804, due to too low thresholds. 

The main reason for the larger error was increased confusion between all 
categories and category 7.2, not site prepared. These errors were small on the 
spectral analysis. 

Using the same spatial post processing that we used in the spectral 
analysis we reduced the error most of the time but not always. After a 4-fill, 
8-fill, 4-shrink, 8-shrink, we eliminated almost all errors in the spectral 
analysis (Section VIII) but with this spectra I -textural analysis (Table XII. 2) 
we increase misidentification error on category 4.1 to 91%, and on category 
2.3 it was about the same (59%) as before post processing. 

The final filling of the image (Table XII. 3) reduced the error rates to 
35% and 31% but did not come close to the 85% classification accuracy of the 
2-band spectral results. This might have been due to the texture function used 
or to the fact that there was little textural distribution between the categories 
in this image. 
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Figure XII. 1 The .82 - .88 micrometer band used for the texture transform. 



Figure XII. 2 The texture transform of Figure XII. 1 with a 3x3 rectangular 

convolution before the texture transform and a 3x3 rectangular 
convolution after. 
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Figure XII. 3 


The texture transform of Figure XII. 1 with no rectangular 
convolution before the texture transform and a 3x3 rectangular 
convolution after. 
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Table XII. 1 The contingency table of the best 2 bond pairs for alpha “ 
beta thresholds of .3 and .021 . 


CONTINGENCY TABLE FOR SAMHA2GDT - 1 SMHA.aABOS - 1 SCALE FACTOR 10»» 
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Table XII. 2 The contingency table of the best 2 band pairs after 4-fill ' 

8-fill, 4-shrink, and 8-shrink operations. ' 
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Table XII. 3 The contingency table of the best 2 band pairs after 4-fill, 
8-fill, 4-shrink, 8-shrink and complete filling operations. 
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XIII Spectral-Textural Analysis; Edit 3 

In addition to the six spectral bands, we provided the feature selector with 
two textural transform bands. The texture bands were created from the .82 - .88 
micrometer spectral bands as before. We used a 3x3 convolution before and after 
textural transform and no convolution before and 3x3 convolution after textural 
transform. The feature selector chose bands .40 - .44 and .588 “ .643 micro" 
meters with .40 - .44 micrometers and no convolution before and 3x3 convolution 
after texture bands for the best 2 band pairs. Figure XIII. 1 shows how the alpha 
and beta thresholds Were chosen in an attempt to minimize the total number of 
reserved decisions and to equalize the number of reserved decisions due to no 
assignment and the number of reserved decisions due to multiple assignments. 

For the best 2 band pairs the alpha threshold was set at .5 and the beta 
threshold at .035. Table XIII. 1 shows the resulting contingency table for the 
best 2 band pairs. There are 49, 130 reserved decisions with 18,083 due to no 
assignment and 31,047 due to multiple assignment. The misidentification error 
rate was 42% and the false identification error rate was 43%. The largest cause 
of error was the misidentification error rate (90%) of category 1.3, shortleaf 
pine, mostly caused by assigning category 1.2, another sublcass of shortleaf 
pine. Post processing with a 4-fill, 8-fill, 4-shrink, 8-shrink and a complete 
filling results in a misidentification error rate of 34% and a false identification 
error rate of 20% (Table XIII. 2). 

The band pairs used for the best 2 along with the .588 - .643 and 
.65 - .69 micrometer band pair were chosen by the feature selector as the 
best 3 band pairs. Figure XIII. 2 shows the graph of the threshold alpha 
against the number of reserved decisions. For the best 3 band pairs the alpha 
threshold was set at .6 and the beta threshold was set at ,042. It is interesting 
to note, that the number of reserved decisions due to no assignment was 25,878, 
and the number of reserved decisions due to more than one assignment was 
26,566 which are very close indicating good thresholds. 

Table XIII. 3 shows the resulting contingency table for the best 3 band 
pairs. The misidentification error rate was 38% and the false identification 
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error rate was 41%. The greatest cause of confusion is the misidentification of 
category 1 .3 and the false identification of category 1 ,3, a subclass of shortleaf 
pine. As with the spectral analysis of edit (Chapter IX), the confusion is 
mostly within class types. Confusion between category 1.2 and category 1.3, 
subclasses of shortleaf pine, and confusion between category 2.4 and category 
2.6 cause most of the error. Figure XIII. 3 shows the best 3 band pairs c|assi“ 
f i cation. 

Post processing with a 4-fili and an 8-fill (Table XIII, 4 and Figure 
XII. 4) did not really change the error rates. The misidentification error rate 
is 40% and the false identification error rate is 44%. A 4-shrink (Table XIII. 5 
and Figure XIII. 5) and an 8-shrink (Table XIII. 6) eliminate almost all of the 
confusion between class types, but the error within class type 1 is still high. 

This confusion within class type 1 was also present in the spectral analysis 
(Chapter IX). The final post processing, a complete filling, resulted in a 
contingency table (Table XIII. 7 and Figure XIII. 6) having a misidentification 
error rate of 25% and a false identification error rate of 29%. 

Note that the results of 4-fill, 8-fill, 4-shrink, 8-shrink, and complete 
f illing operations for the best 3 band pairs (Table XIII, 7 and Figure XIII. 6) 
and the results of 4-fill, 8-fill, 4-shrink, 8-shrink and complete filling 
operations for best 3 band pairs using the spectral analysis (Chapter IX, Table 
IX .7 and Figure IX .8) shows less error in the spectral results. Yet, comparison 
of Figure XIII. 6 and Figure IX .8 show that the figure from the texture analysis 
is actually truer to the timber stand and compartment map for edit ^3 than the 
spectral figure. It seems this is due to the area covered by the ground truth 
overlay (Figure IX. 2), so that more ground truth would have resulted in better 
classification accuracy for the texture analysis. 
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Figure XIII. 1 Number of reserved decisions as a funcfion of 

probability threshold alpha for best 2 band pairs 
with texture for edit ^3 
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Table XIII. 1 The contingency table of the best 2 band pairs for alpha “ 
beta thresholds of .5 and .035. 
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Table XIII. 2 The contingency table of the best 2 band pairs after 4-fill, 
8-fill, 4-shrink, 8-shrink and complete filling operations. 
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Figure XIII. 2 Number of reserved decisions as a function of 

probabil ity threshold alpha for best 3 band pairs 
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Table XI11.3 The conHngency table of the best 3 band pairs for alpha 
beta thresholds of .6 and .04-2. 
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Table XIII.4 The contingency table of the best 3 band pairs after 4-fill 
and 8-fill operations. 
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Figure XIII. 3 


The classification of the three best band pairs for alpha 
beta thresholds of .6 and .042. 



Figure XIII. 4 The classified image of Figure XIII. 3 after 4-fill and 8-fill 
operations. 
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Table XIII. 5 The contingency table of the best 3 band pairs after 4“fill, 
8-fiil, and 4^shrink operations. 
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Table XIII. 6 The contingency table of the best 3 band pairs after 4“f|ll, 

8-fill, 4-shrink and 8-shrink operations. 
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Table XIII, 7 The contingency table of the best 3 band pairs after 4“fill, 
8"fill, 4-shrink, 8~shrink and complete filling operations. 



Figure XIII. 5 




Figure XIII. 6 


The classified image of Figure XIII. 3 after 4-fill, 8-fill 
and 4-shrink operations. 



The classified image of Figure XIII. 3 after 4-fill, 8-fill, 
4-shrink, 8-shrink, and complete filling operations. 
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XIV Conclusions 

Use of textural features and spatial post processing has been shown to cut the 
average classification error to less than half its initial value while tending to increase 
and equalize the equally weighted average misidentification error and equally weighted 
false identification error. The classified images resulting from spatial and textural pro” 
cessing have a more cartographic map-like quality than the typically salt and pepper 
classified images using no textural features or spatial post-processing. 

The simultaneous decrease in average classification error and increase in both 
equally weighted average misidentification error and false identification error means 
that more pixels whose true category identification is of a frequently occurring cate- 
gory get reassigned correctly by the spatial post-processing than of an infrequently 
occurring category. This is a natural consequence of the fact that the spatial pro- 
cessing is more of a syntactic operation than a semantic one. Spatial processing 
operations which use category labels intead of just sameness or difference of category 
labels could be designed which do not favor the larger categories over the smaller 
categories. 

Because of the strong interaction between average error, average misidenti- 
fication error, average false identification error, and classified image appearance 
and complexity, it is clear that further work can bear much fruit by analysis of these 
interaction effects. In particular, we recommend that textural and spatial post- 
processing concepts be developed using classified image's local neighborhood contexts 
as the independent variable and classification error as the dependent variable. 
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APPENDIX 1 

Textural Transform Programs 



PROGRAM DESCRIPTION 


Table of Contents 

I . User Interaction 
II. Internal Program Description 
111, Non-Standard Subroutines 
IV. Subroutine Documentations 
V. Listing 
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User Interaction 


User parameters are input by the routine TXINPT which asks for para- 
nieLers by name: 

NFUNC = 1 use subroutine FUNC 1 

NFUNC = 2 use subroutine FUNC 2 

NFUNC = 3 use subroutine FUNC 3 

NFUNC = A use subroutine FUNC k 

NFUNC = 5 use subroutine FUNC 5 

NDIS = distance between spaces of neighboring cells 

IBOUT = logical unit number to output error messages 

PCLCT = percent of point to count In FPLXIT 

FILNMP = input filename 

FILNMQ, = output filename 



I nternal Program Descript ion 


The texture programs are set up so that after the call to TXINPT at 
the beginning of execution the user does not interact anymore with the 
computer. 

The user must know which function he wants to use. The user inputs 
1, 2 , 3, Of 5 corresponding to FUNCl, FUNC2, FUNC3, FUNC4, and FUNC5 
respectively; where: 

FUNCl - computes the sum probability feature of the image 

FUNC2 - computes the gradient entropy feature of the image 

FUNC3 ~ computes the entropy feature of the image 

FUNC4 - computes the gradient feature of the image 

FUNC5 “ prepares normalized lex arrays which have been equal proba- 
bility quantized according to their diagonal elements. 

The parameter NDIS is the distance between spaces of neighboring cells. 

The texture transform works on the co-occurrence of grey levels on neigh- 
boring cells. Each cell has a 0° neighbor, a 90° neighbor, a 135° neighbor, 
and a 45° neighbor. This covers all the cells, since a cell's 180° neighbor 
has that cell as a 0° neighbor. Thus, for each grey level, there is a count 
of the co-occurrences of grey levels as one of the four specified neighbors. 

The parameter NDIS is the distance the algorithm gives to look for the 
neighboring cells. If the user wants to perform the texture transform using 
all co-occurrence counts, then the parameter PCLCT should be 1,00. If the user 
only wants to count 80 % of the cells, the PCLCT should be set to 0.80, and 
so on. 

The mainline TXJDM calls the ASCII i/0 routine TXINPT for input para- 
meters. The TXJDM transfers control to TXTMN. This routine sets up the work 
area, allocating core to those arrays that need it. FPLX1T is then called 
to compute the lex arrays where: 

LEXi - array containing count over all grey levels of vertically 
adjacent (90-degree) neighbor; 

LEX2 - array containing count over all grey levels of horizontally 
adjacent (O-degree) neighbor; 

LEX3 - array containing count over all grey levels of left diagonally 
adjacent (135-degree) neighbor; 

LEX4 - array containing count over all grey levels of right diagonally 
adjacent (45-degree neighbor. 



When these counts have finished, control is returned to TXTMN, which 
transfers control to the appropriate function as specified by the user. 

The FUNC array which is passed as an arqurnent to the FUNC routines is 
equivalanced to the lex arrays. For example: 

FUNC (1,1) = LEXl (1 ) 

FUNC (1,2) = LEX2 (l) 

FUNC (1 ,3) - LEX3 (1 ) 

FUNC (1,4) = LEX4 (1 ) 

After the appropriate function has been applied, control is again returned 
to TXTMN. TXTMN then calls in PLXIT. PLXIT reads in the image data and 
determines the corresponding eight neighbors and applies the texture transform. 

Let x be the set of resolution cells of an image I (by row- 
column coordinates). Let G be the set of grey tones possible to appear on 
image i. Then I: Z x Z ->■ G. Let R be a binary relation on Z x Z 
pairing together all those resolution cells in the desired spatial relation. 

The co-occurrence matrix P, P: G x G [0,1], for image 1 and binary 

relation R is defined by 

P(i,j) = 

The textural transform J,J: Z^ x Z^ (-“,“), of image I relative to function » 

f, is defined by 

J(y>x) = ,, J V 2 f [P( i (y,x) , I (a,b) ) ] 

(a,b)cR(y,x) 

Assuming f to be the identity function, the meaning of J(y,x) is as 

follov^;5. The set R(y,x) is the set of all those resolution cells in 

Z X Z in the desired spatial relation to resolution cell (y,x). For 
r e 

any resolution cell (a , b) cR (y , x) , P ( I (y ,x) , I (a , b) ) is the relative frequency 
by which the grey tone l(y,x), appearing at resolution cell (y,x), and the 
grey tone l(a,b), appearing at resolution cell (a,b), co-occur together in 
the des i red spati al relation on the entire image. The sum 

P(l(y,x), l(a,b)) 

(a,,b)r.R(y , x) 

is just the sum of the relative frequencies of grey tone co-occurrence over 


l4l 


IS 



{((a,b),(c,d))eR|l(a,b) = i and I (c,d) = J } 


all resolution cells in the specified relation to resolution cell 


(y,x). 


The factor 


1 


the reciprocal of the number of resolution cells in the 


desired spatial relation to (y,x) is just a normalizing factor. 

These data values are then written out in the corresponding place on 
the output texture transformed image. When PLXIT exits, the texture trans- 
form has been created. Control goes to TXTMN, which exits back to the 
mainline TXTDM. This program returns to the beginning and brings back the 
ASCII I/O routines to get the parameters for the next texture transform. If 
none are desired, a carriage return will terminate the processing. 

All ASCII I/O on our PDP-15 is 5/7 ASCII in double integer words. The 
PDP-15 has 36 bits in one double integer word. 

See Figure 1 for the program flow. 
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III. Non-Standard Subrout i nes 
ADJl 

ADJ2 Dynamic core allocation routines 
ADJ3 

The program can allocate memory by performing what essentially amounts 
to a dynamic Fortran equivalence and dimension 
KDPUSH - ignore (delete) 

Error stack processing used In KANDI DATES. 

SDKINL - KANDIDATS sequential file opener 

Opens files for KANDIDATS routines. Uses Seek and Enter (STANDARD 
Fortran routines) and can be modified to fit your file structure. 

SKPDSC - skip descriptor records 

KANDIDATS creates descriptor records, containing processing history 
information, before the image date. Since the file is sequential, these 
must be skipped. If the user has random access on images, this can be 
ignored. If not, be sure that image record numbers are advanced to first 
image data record. 

IMTRXP - Matrix print-out routine 

Any standard Matrix print routine will work. 

SREAD - sequential read (uses Fortran reads) 

SWRITE - sequential write (uses Fortran write) 

Starting on the next page is an explanation of the "ADJ routines 
and several ideas on how to get around them. 
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The included program segment may be compared to the following example: 

1. INTEGER ARRAY (500), X(l) 

2. REALY(l) 

3. INTEGER XSIZE, YSIZE, YSTART, TTLSZE 

4. READ (5, 100) XSIZE, YSIZE 

5. TTLSZE = XSIZE + YSIZE 

6. IF (TTLSZE .GT. 500) CALL ERROR 

7. YSTART = XSIZE + 1 

8. CALL ADJl (X, ARRAY (1)) 

9. CALL ADJl (Y, ARRAY (YSTART)) 

10. ...TASK CODE 

11. STOP 

12. END 

Within the task code, X and Y may be referenced as vectors with respective 
types, Integer and real. In addition references to X will access the first 
XSIZE elements of ARRAY and references to Y will access the last YSIZE 
elements of ARRAY. 



If X and Y are used only in contexts that functions may be used in, then the program 
segment may be recoded using statement functions* (Check your particular imple- 
mentation of FORTfl^.N for applicability.) 

1. INTEGER ARRAY (500) 

2. REAL RL 

3. INTEGER XSI2E, YSIZE, TTLSZE 

4. X(I) = ARRAY (I) 

5. Y(I) = RL(ARRAY (1 + XSIZE)) 

6. READ (5, 100) XSIZE, YSIZE 

7. TTLSZE XSIZE + YSIZE 

8. IF (TTLSZE .GT. 500) CALL ERROR 

9. ... TASK CODE 

10. STOP 

n. PND 

Where the function RL is coded as follows. 

1. REAL FUNCTION RL(ARG) 

2. REAL ARC 

3. RL = ARG 

4. RETURN 

5. END 

V/ithin the task code, X and Y will have respective types integer and real 
and will access those specified locations of ARRAY. 

However, X and Y may only be used as functions. 


1^5 


. In the context of subroutine colls, adjustable dimensions is a standard feature 
of FORTRAN as In the following example; 

1. INTEGER ARRAY (500) 

2. INTEGER XSIZE, YSIZE, YSTART, TTLSZE 

3. READ (5, 100) XSIZE, YSIZE 

4. TTLSZE = XSIZE + YSIZE 

5 . IF (TTLSZE .GT. 500) CALL ERROR 

6. YSTART = XSIZE + 1 

7. CALL SUB (ARRAY^l), ARRAY (YSTART), XSIZE, YSIZE) 

8. STOP 

9. END 

Where SUB is coded os follows: 

1. SUBROUTINE SUB (X, Y, XSIZE, YsizE) 

2. INTEGER XSIZE, YSIZE 

3. INTEGER X(XSIZE) 

4. REAL Y(YSIZE) 

#• 

5. ... TASK CODE r. 

6 . RETURN 

This approach necessitates a division of storage allocation code and task code. 
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Alternatively X and Y may be dimensioned independently and given a 
reasonable but sufficient size. 


1 . INTEGER X (250) 

2 . REAL Y(250) 

3. READ (5, 100) XSIZE, YSIZE 

4. IF (XSIZE .GT. 250). OR. 

5 . (YSIZE .GT. 250) CALL ERROR 

6. ... TASK CODE 

7. STOP 

8. END 




Check the output of the FORTRAN compiler being used. 

If the compiler generates and uses dope vectors it v/ould be possible to produce 
user written ADJ routines. 

Keep in mind that all recoding must preserve the size, shape, type and 
usage of the involved data elements. 


V. Subroutine Documentations 


ORIGINAL PAGE IS 
OF POOR QUAUTYi 
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GENERAL MATRIX PRINTOUT PROGRAM 


PROGRAM TITLE; 

DATE OF LISTING: 
PROGRAMMER: 
DOCUMENTED BY: 
PROGRAM LANGUAGE; 
COMPUTER REQUIRED: 
PURPOSE: 


CALLING SEQUENCE; 


INPUT ARGUMENTS: 


SUBROUTINE IMTRXP 
February 13, 1973 
Dinesh Goel 
Dinesh Goei 
FORTRAN IV 
PDP 15/20 

Th is subroutine divides an integer matrix into sections 
suitable for printer output and prints the matrix with 
matrix title, column designation, row designation, and 
column and row labels. 


CALL IMTRXP (lA, NROW, NCOL, NRWDIM, TTLl, TTL2, 
TTL3, CLBL, RLBL, ISTR) 


lA 

NROW 

NCOL 

NRWDIM 

TTLl 

TTL2 

TTL3 

CLBL 

RLBL 

ISTR 


Input array of matrix to be printed out. 
Number of rows in the printed matrix. 

Number of columns in the printed matrix. 

Row dimension of the entire matrix which is 
stored by columns. 

Matrix title of 13 words. 

Column title of 2 words. 

Row title of 2 words. 

Array of column labels. 

Array of row lables. 

This is an option 

if 1, matrix will be printed as such, 
if 2, transposed matrix will be printed out. 
if 3, matrix is assumed to be symmetric having 
long to short storage in !A. 
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if 4, matrix Is assumed to be summetric having 
short to long storage in lA. 


OUTPUT ARGUMENTS: 

None . 

OTHER PARAMETERS AND ARRAYS: 

IROW Array for any one row of matrix as finally 
printed out. 


COMMENTS: 

If the printed matrix has large number of columns which can 
not fit on one page of printer output, it will be separated 
into blocks, each of which is small enough to fit on one 
page. The rows are printed in the blocks of 5. This 
program takes only the integer numbers, for real numbers 
RMTRXP can be used. File code 17octal has been used 
for printing the matrix which must be assigned to teletype 
or IBM printer in the beginning as desired. 
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Sequential Read 


PROGRAM TITLE; SREAD 

VERSION : B 


DATE: 

UPDATE: 

AUTHOR; 

DOCUMENTED BY: 
PROGRAM LANGUAGE; 
IMPLEMENTED ON; 
PURPOSE; 


June 22, 1973 
April 29, 1975 
Robert M. Haro lick 
Robert M. Haro lick 
FORTRAN IV 
PDP 15/20 


This subroutine reads a set of lines on a file of single image 
data stored in standard bit compacted form, SREAD assumes 
that SDKINL has already been called to open the file on IDAT. 


ENTRY POINT; 

SPREAD (IDAT, lARRAY, IDY, IDENT, lEV, ERRRET) 


ARGUMENT LIST; 

IDAT 

lARRAY 

IDY 


IDENT 

lEV 


the file code on which the image resides 

2 ■'dimensional array (row x column) which 

subimage is returned in 

number of records in subimage. The number 

of rows and columns for a record will be taken 

from IDENT (14) and IDENT (13), respectively 

identification array of 20 words. 

integer event variable 

lEV = 1 success 

lEV = “2001 illegal file code 

lEV = “2006/1 0/too small 

lEV = “2007 EOF 

lEV = “2009 READ ERROR 

lEV = 2012 illegal data mode 
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ERRRET alternate return taken if an error occurs 

SUBROUTING REQURED; 

ADJl 

KDPOP 

KDPUSH 

UNPACK 

COMMENTS: 

lARRAY must be a two dimensional array 

IDENT (13)* IDENT (14) must not be greater than 256 
unless the user has a block data program to allocate more 
memory to labeled common area /\ 0 / , 
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Skip Descripfor Records 


PROGRAM TITLE; 
VERSION: 

DATE: 

UPDATE: 

AUTHOR: 

DOCUMENTED BY: 
PROGRAM LANGUAGE: 
IMPLEMENTED ON: 
PURPOSE: 


SKPDSC 

A 

July 10, 1973 
October 15, 1974 
Robert M. Horolick 
Robert M, Horolick 
FORTRAN IV 
PDP 15/20 


This progrom skips the descriptor records of imoges stored in 
stondord bit compocted form. SKPDSC ossumes thot SDKINL 
hos been colled previously* 


ENTRY POINT: 

SKPDSC aOATP, IDENT, lEV, ERRRET) 


ARGUMENT LIST: 

IDATP 

IDENT 

lEV 


ERRRET 

SUBROUTINES REQUIRED: 
KDPOP 
KDPU5H 


file code on which the image resides, 
identification array of 20 words for the image. 
= 1 success 

= -2001 illegal file code 

= -2007 EOF 
= “2009 read error 

Alternate return taken if error occurs 
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Sequential Disc Initializer 


PROGRAM TITLE: SDKINL 

VERSION: B 


DATE; 

UPDATE: 

AUTHOR; 

DOCUMENTED BY; 
PROGRAM LANGUAGE: 
IMPLEMENTED ON: 
PURPOSE: 


June 30, 1973 
October 15, 1974 
Robert M. Haralick 
Robert M. Haralick 
FORTRAN IV 
PDP 15/20 


This subroutine initializes a PDP 15/20 sequential disc file 
for input or output. The file is used to store image data in 
standard bit compacted form. The number of data words will 
be written in a logical record of at least 20 and the number 
of bits per data word should not be more than 18. 


ENTRY POINT: 

' SDKHNL: (iDAT, FiLNM, iuENT, iRDVVRT, iEV, ERRRET) 


ARGUMENT LIST: 

IDAT 

FILNM 

IDENT 

IRDWRT 


lEV 


file code on which file resides, 
array containing the fil name, 
identification array of 20 words, 
read/write Indicator. 

IRDWRT =1 initialize as input file. 

IRDWRT =2 initialize as output file 
integer event variable. 

= 1 Success 

= -2001 Illegal file code 
= -2002 Number of bits per point has a 
illegal value 

= -2003 Frame coordinate and image dimension 
infonnation not specified in-Ident-array 
= -2004 Illegal request 
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ERRRET 

SUBROUTINES REQUIRED: 
ENTER 
FSTAT 
ICEIL 
KDPOP 
KDPUSH 
MAXj2f 
SEEK 


= -2005 file does not exist 

= “2011 illegal min/max/lNjZL/nbits combination 
= "2012 illegal data mode 

ALTERNATE RETURN TAKEN IF A ERROR OCCURS 
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PROGRAM TITLE: 
VERSION: 

DATE: 

UPDATE: 

AUTHOR: 

DOCUMENTED BY; 
IMPLEMENTED ON 
PURPOSE: 


ENTRY POINT: 


ARGUMENT LIST: 


Sequential Write 

SWRITE 

C 

June 22, 1973 
October 15, 1974 
Robert M. Horalick 
Robert M. Horalick 
: PDP 15/20 

This program writes a set of lines or a file of single image data stored 
in standard bit compacted format. SWRITE assumes that SDKINL 
has also been called to initialize the file on IDAT'. 


SWRITE 0DAT,1ARRAY, IDY, IDENT, lEV, ERRRET) 


IDAT 

lARRAY 

IDY 


IDENT 

lEV 


ERRRET 


file code on which 1 image resides. 

2-dimensional array (row x column) 

in which subimage is transferred to program. 

number of records for subimage. The number 

of rows and columns for a record will be taken 

from IDENT (13) and IDENT (U). 

identification array of 20 words. 

integer event variable 

lEV = 1 success 

1EV = -2001 illegal file code 

lEV “ “2006 /lO/ too small 

1EV = -2007 EOF 

1EV = -2008 WRITE ERROR 

lEV = -2012 illegal data mode 

ATTENATE RETURN TAKEN IF ERROR OCCURS 
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SUBROUTINE REQUIRED: 
ADJl 


KPPUSH 
KDPOP PACK 

COMMENTS: 

lARRAY must be a two dimensional array. 

IDENT (13)* IDENT (14) must not be greater than 256 
unless the user has a block data program to allocate more 
memory to labeled common area /lO/. 
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OEIGmAi; PAGE IS 
OE! POOR QUAUTy 
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CERROR 

C 

c 

c 


E-R-R-O-R 


ASCII ERROR I/O FOR TEXTURE PROGRAM 


C 

c 

c 

c 

c 

c 

c 

c 


IJ 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


u 

c 

w 405 
304 



311 

400 

200 


PROGRAM TITLE 

VERSION 

AUTHOR 

DATE 

UPDATE 

PROGRAM LANGUAGE 
IMPLEMENTED ON 
DOCUMENTED BY 


ERROR 

A 

CH IN-HUANG CHEN 
FEBRUARY 1975 

FORTRAN IV 
PDP 15 

CH IN-HUANG CHEN 


PURPOSE 


THIS ROUTINE. TELLS THE USER EITHER LSTID OR TXTMN I 
IN ERROR ON . DAT SLOT IOU( OR IBOUT 


ENTRY POINT 
ARGUMENT LIST 
I ERR 

lEV 
I OUT 
I BOUT 


ERRORdERR. lEV, I0UT3 iBOUT) 

PARAMETER USED TO DETERMINE EITHER LSTID 
OR TXTMN IS IN ERROR 
INTEGER EVENT VARIABLE 
ERROR MESSAGE OUTPUT . DAT SLOT 
ALTERNATE ERROR MESSAGE OUTPUT . DAT SLOT 


SUBROUTINE ERROR (I ERR. lEV. I OUT. IBOUT) 
DOUBLE INTEGER FDATE<3) 


GU FU 'l3U4’. 3lU). I ERR 
CALL ABATE (FDATE) 

WRITE (I OUT. 405) FDATE 
FORMAT/ IX. 3A5) 

I F ( I BOUT. NE. I OUT ) WR I TE C I BOUT . 405 ) FDATE 
GO TO 200 

WRITEC lOUT. 305) lEV 

I F (I OUT. NE. I BOUT ) WR I TE ( I BOUT . 305 ) I E V 
FORMAT(-' LSTID ERROR-’.IS) 

GO TO 400 

WR I TE < I BOUT . 3 1 1 ) I EV 

I F ( I BOUT. NE. I OUT ) WR 1 TE ( IBOUT .311) I EV 
FORMAT (••• TXTMN ERROR IEV==-L. 15) 

GALL CLOSE (I EiOUT) 

RETURN 

END 
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F-P-L-X-I-T 


CFPL 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


u 


XIT 


PRiJijiRAl’l 

TITLE 


SUEROUT I ME FPLX I T 


PROGRAMMER A. SINGH NOVEMBER 1.972 

UPDATE ROBERT M HARALICK FEBRUARY 1974 

GE MONAGHAN 9/20./74 

RM HARALICK 10./ 10/74 

CH IN-HUANG CHEN 2/22/75 

PURPOSE ADD PCLCT IN ARGUMENT LIST 

CHANGE LEX ARRAY TO SINGLE INTEGER 
ADD OVERFLOW CHECK ON LEX ARRAY 


DOCUMEN- 

TATION 


A. SINGH 


COMPUTER 

REQUIRED 


ANY 


i; 


PROGRAM 

LANGUAGE 


FORTRAN IV 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PURPOSE 


METHOD 


CALLING 

SEQUENCE 

ARGUMENTS 


FPLX IT COMPUTES THE FOUR NEIGHBOUR GRAY TONE 
MATRICES LEXl, LEX2.. LEX3 AND LEX4 FOR ANGLES 90,0,135 
AND 45 DEGREES RESPECTIVELY. IT WORKS FOR ALL 
DISTANCES, 

FPLX IT CHECKS THE GRAY LEVELS OF THE NEIGHBOURS OF 
A CELL, AND INCREMENTS THE CORRESPONDING ELEMENT IN 
THE ASSOCIATED LEX ARRAY. THE NEIGHBOURS UNDER 
CONSIDERATION ARE A DISTANCE ?<D?< AWAY, WHERE ?<DS< 

IS THE DISTANCE FOR THAT RIJN OF FPLX IT. 

CALL FPLXIT( IDATI, IDATA, LEX 1 , LEX2, LEX3, LEX4, IPT, IDENT, 
MMl, PCLCT, lEV, lERRl) 


U 

C 

C 

C 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


IDATI INPUT FILE CODE 

IDATA SCRATCH ARRAY FOR MMl LINES OF THE IMAGE 


LEX 1 , LEX 

IPT 

IDENT 

MMl 

PCLCT 

lEV 


lERRl 


IMAGE. 

LEX3 AND LEX 4 ARE THE FOUR LEX ARRAYS 
FOR THE GRAY TONE MATRICES. 

ARRAY WHICH CONTAINS THE POINTERS FOR 
THE IDATA ARRAY. 

IDENTIFICATION ARRAY FOR THE IMAGE 
SPATIAL DISTANCE -H 1 
PERCENT OF LINES COUNTED 
INTEGER EVENT VARIABLE 

IEV=-5011 IF NLIMPPL OR NUMLIN IS LESS 
THAN TWICE SPATIAL DISTANCE PARAMETER. 
IEV=-5010 IF LEX ARRAY IS OVERFLOW 
ALTERNATE RETURN TAKEN IF ERROR OCCURS 
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c 

PARAMETERS 

NUMLIN 

NUMBER OF LIMES IN THE IMAGE 


c 


NUMPPL 

NUMBER OF POINTS PER LINE IN THE IMAGE 


c 


I MAX 

LARGEST GRAY TONE 

- 

c 


IMIN 

LEAST GRAY TONE 


c 


LEAST 1 

=IMIN-1. LEAST 1 IS USED FOR NORMALISING 


c 



THE GRAY TONES. 

- 

c 


NOBL 

NUMBER OF GRAY TONES 


c 



N0BL=IMAX-IMIN+1 


c 


NBUBL 

SIZE OF A LEX ARRAY 

- 

c 



NBUBL=NOBL^«- ( NOBL+ 1 ) /2 


c 

)-■ 

INPUT AND 

IMAGE READ IN FROM FILE (02). 

- 

c 

OUTPUT 




c 

c 

RETURNS 

NORMAL 

AND ALTERNATE ERROR RETURNS 


c 

SUBPROGRAMS 

INDEX 



c 

REQUIRED 




c 

CALLED BY 

TXTMN 



c 

c 

COMMENTS 

FPLX IT 

WORKS FOR ALL SPATIAL DISTANCES. IT DOES THIS 


c 


BY HAVING NDIS+1 LINES OF I DATA IN CORE. WHERE NDIS 


c 


IS THE 

SPATIAL DISTANCE PARAMETER. 



SLIBR0IJTI^'4E FPLX I T ( IDATI , IDATA, LEXl. LEX2, LEX3. LEX4, IPT; 
IIDEMT. Mill, IMGNO, PCLCT.. I ARRAY. lEV. lERRl ) 

DOUBLE INTEGER INT/ LEX i. LEX2. LEX3. LEX4 
DIMENSION IDATA< I, 1 ). LEXl ( 1) . LEX2( 1). LEX3( 1) . LEX4( 1 ) . IPT( 1 ) 
D I MENS I ON I DENT { 2.0 ) . I ARRAY (1.1,1) 


C 


U 

G 


I DATA ( NUMPPL, MM 1) . I PT ( MM 1 ) 

STACK SUBROUTINE NAME IN ERROR STACK 
CALL KDPUSH( •-FPLXl--. •-T-' ) 


SET PARAMETERS 

NUMPPL=IDENT(6) 

NUMLIN=IDENT(7) 

IMIN=rDENT( 15) 

IMAX=IDENT( 16) 

LEAST1=IMIM-1 
NOBL= I MAX --LEAST 1 
NBUBL=NOBL-i«- ( NOEL-h 1 ) /2 



DO 14 1 = 1. NBLIBL 


INITIALISE THE LEX ARRAYS TO ZERO 
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LEX 1(1) =0 
LEX2( I )=0 
LEX3C I 

14 LEX4(I)=0 


CHECK IF SIZE GF IMAGE IS TOO SMALL. 
RELATIVE TO THE DISTANCE PARAMETER MM 

MM=MM1-1 

MM;?-MM-!<-2 

XEV-==-5011 

IF(MUMPPL. LT. MM2. OR. NUMLIN. LT. MM2) GO TO 9999 

NUMPMM=NUMPPL-MM 

MLIMLMM-NUMLIN-MM 

READ IN THE FIRST MMl LINES OF THE IMAGE 
AND SET UP POINTERS 

DO 110 IY=1.MM1 
IPT(IY)=IY 

CALL RREADdDATI. lARRAY. IMGNO, lY. 1. IDENT. lEV. ERRl ) 

DO ill LY=1,NUMPPL 
111 IDATA(LY/ IY) = IARRAY( 1,. 1, LY) 

110 CONTINUE 


SETTING UP POINTERS FOR THE FIRST AND 
LAST ROWS OF THE IMAGE ARRAYS 

I ST- I PT Cl) 

LST=IPT(MM1) 

GO THROUGH ALL BUT MM ROWS OF IMAGE 

NOVFLO=131017 

INT=3S56347S31 

NEXT-MMl+1 

DO 105 LCNT = NEXT. NUMLIN 
IFCRCMC INT). GT. PCLCT) GO TO 105 

SKIP LINES RANDOMLY BY USING RANDOM NUMBER 
GENERATOR RCM EXTERNAL FUNCTION 


GO THROUGH EACH ROW MM TIMES. THE FIRST 
SET OF MM COLUMNS ARE HANDLED SEPARATELY 

DO 120 IRW==1,MM 

OF pooB QUABrrx 
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IRM=IRW+MM 


SET I, L; J AND K EQUAL TO THE 
(NORMALISED) VALUES OF GRAY TONES OF 
RESOLUTION CELLS IN POSITIONS A. B, C 
AND D AS IN THE DIAGRAM ~ 

A C 
B D 

WHERE A INITIALLY IS THE UPPER LEFT 
CORNER CELL. THE CELLS ARE A DISTANCE 
MM APART. 


I = I DAT A ( I RW .. I ST ) -LEAST 1 
L= I DATA ( I RW , LST ) -LEAST 1 
K= I DATA ( IRM ■ LST ) -LEAST 1 
U= I DATA ( I RM. I ST ) -LEAST 1 

PUT THE TWO DIMENSIONAL INFORMATION 
INTO ONE DIMENSIONAL FORM. THE FUNCTION 
NEEDED TO CONVERT A DOUBLE SUBSCRI-TED 
ARRAY. IMM(X.Y).. INTO A SINGLE 
SUBSCRIPTED ARRAY. IMM(Z). IS OF THE 
FORM G ( X ) + F ( Y ) . WHERE G ( X ) = ( X - 1) -«-X /2 
AND FCY) = Y. THEREFORE 
Z = (X-l)*X/2+Y 

THIS IS DONE IN THE PROGRAM BY THE 
EXTERNAL FUNCT I ON I NDEX ( X . Y ) . 

SINCE THE ORDER OF OCCURRENCE OF THE 
GRAY TONES BELONGING TO A RESOLUTION 
CELL PAIR IS IMMATERIAL. THE ARRAYS ARE 
SYMMETRIC. WE LET THE LARGER OF THE TWO 
HAVE THE FIRST SUBSCRIPT. I. E. . THE ARRAY 
IS STORED IN LOWER TRIANGULAR FORM. THE 
ORDER OF THE SUBSCRIPTING IS AS FOLLOWS - 
IMM(l.I) = IM(1). 

I MM (2.1) = IMM(2). 

I MM (2. 2) =•■ IMMC3). 

IMM(3. 1 ) = IMM(4) . 


I MM ( NOBL. NOBL ) = I MM ( NBUBL ) . 

THE SCANNING PROCEDURE. THAT IS THE 
METHOD BY WHICH THE PAIRWISE COMPARISONS 
ARE MADE. IS DESCRIBED BELOW FOR THE 
GENERAL CASE. 

CONSIDER A RESOLUTION CELL WITH SPATIAL 
COORD I NATES ( M . N ) . AND CALL TH I S CELL I . 
THE SCANNING OPERATION BEGINS IN THE 
UPPER LEFT HAND CORNER OF THE IMAGE AND 
IT THEN PROCEEDS BY COMPARING THE GRAY 
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TONE OF WITH AT MOST FOUR GRAY TONES 

OF ITS MEIGHEOLIRING RESGLUriOW CELLS. 
THAT NEVER NEEDS TO CONSIDER MORE 

THAN FOUR NEIGHBOURS CAN BE SEEN FROM 
THE DIAGRAM OF THE SEARCH PATTERN SHOWN 
BELOW — 


I J 
M L. K 

ON A GIVEN ITERATION, ?/I?c WILL LOOK FIRST 
AT ITS VERTICAL NEIGHBOUR > NEXT 

AT ITS HORIZONTAL NEIGHBOUR THIRD 

AT ITS LOWER RIGHT NEIGHBOUR (S<K?4<) AND 
FOURTH AT ITS LOWER LEFT DIAGONAL 
NEIGHBOUR (S^M?/.). Zf.U< THEN MOVES INTO 
THE POSITION OF THE RIGHT RESOLUTION 
CELL OF THE PREVIOUSLY SCANNED FIRST 
ROW (THE POSITION OCCUPIED BY 
THE OPERATION IS REPEATED UNTIL ALL 
NEIGHBOURING PAIRS OF RESOLUTION CELLS 
HAVE BEEN EXAMINED. THE PROCEDURE IS 
FURTHER REPEATED FOR CELLS SKIPPED OVER 
IF THE SPATIAL DISTANCE IS GREATER THAN 
ONE, TILL ALL CELLS HAVE BEEN EXHAUSTED. 


T1..,'-INDEX( T, L) 

COUNT VERTICALLY ADJACENT (90-DEGREE) 
ME I GHBOLIR 

LEX.1. ( IL)=LEX1 ( IL ) + l 
I J=TMDEX( I, J) 

COUNT HORIZONTALLY ADJACENT (0-DE+REE) 
NEIGHBOUR 

LEX2( TJ)=LEX2( IJ)+1 
Ik-IMDEX(I.1<) 

COUNT LEFT DIAGONALLY ADJACENT 
(13S-DEGREE) NEIGHBOUR 

I EX3aK)“LEX3aK)+ l. 

NOW ITERATE DOWN THE ROW 

DO ISO N-IRM, NUMPMM, MM 

WMM=N+MM 

I=J 
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M=L 

L=K 

C 

J=IDA7A<NMM; I ST) -LEAST 1 
K= I DATA ( NMM, LST ) -LEAST 1 

0 

IL=--INDEX( I. L) 

C 

C COUNT VERTICALLY ADJACENT (90-DEGREE) 

C NEIGHBOUR 


U 


c 

C- 


LEXl ( IL)=LEX1 ( ID+l 


IJ=INDEX( I, J) 


C 

C 

C 

C 

C 

iC 

c 

c 

c 

c 


:C 



LEX2 ( I J ) =LEX2 < IJ ) + 1 
IK=INDEX( I, K) 


LEX3< IK)=LEX3( IK)+1 
IM=INDEX( I, M) 

LEX4( IM)=LEX4( IM)+1 
130 CONTINUE 


j 

i C 

I c 
! C 
i C 

i c 


; c 
I c 
I c 
’ c 


I=J 

M=L 

L=K 

IL=INDEX( I, L) 


LEXl < IL)=LEX1 ( ID +l 


I M= INDEX ( I, M) 


COUNT HORIZONTALLY ADJACENT (0-DE+REE) 
NEIGHBOUR 


COUNT LEFT DIAGONALLY ADJACENT 
( 1 35-DEGREE ) NE I GHBOUR 


COUNT RIGHT DIAGONALLY ADJACENT 
( 45-DEGREE ) NE I GHBOUR 


COMPUTE THE LAST SET OF MM COLUMNS 
SEPARATELY 


COUNT VERTICALLY ADJACENT (90-DEGREE) 
NEIGHBOUR 


COUNT RIGHT DIAGONALLY ADJACENT 
(45-DEGREE) NEIGHBOUR 


166 



LEX 4 (IN) =LEX 4 < I M 3 + 1 


C 

120 COMFINUE 
C 


SHIFT THE F'O INTERS FOR THE TWO ARRAYS. 
THIS IS DONE BY A CYCLIC ROTATION. 

THE POINTER ARRAY IPT IS SUCH THAT AT ANY 
TIME THE ITH LOCATION OF IPT CONTAINS 
THE POINTER TO THE iTH POSITION OF THE 
LINE IN I DATA OR JDATA ARRAY. FOR 
EXAMPLE.. IF IPT (2) =4 THEN THE. FOURTH LINE 
OF THE PHYSICAL JDATA ARRAY IS ACTUALLY 
THE SECOND LINE. AT THAT MOMENT. 


IF(LCNT, EQ. NUMLIN) GO TO 105 

ROTATE IN A CYCLIC MANNER 

ITEMP=IPT( 1 3 
DO 135 IB=:1.MM 
135 IPT(IB3=-IPT( IB+1 3 
IPT ( MM 1 3=1 TEMP 

SET UP THE POINTERS TO THE FIRST AND 
L.AST ROWS OF THE TWO IMAGE ARRAYS 

IST=IPT(13 


lst=ipT':mmi) 

c 

C READ IN A NEW LINE INTO THE I DATA ARRAY 


112 


105 

106 


IJ 

C 

C 

c 


c 


c 


c 


CALL RREAD( IDATI.. I ARRAY. I MONO. LCMT. 1. I DENT 
DO 112 LY=1.NUMPPL 
I DATA L Y. LST 3 = I ARRAY i 1 .. 1 .. L Y 3 
IF(LEX1( IL3. GT. NOVFLO) GO TO 106 
IF(LEX2aJ3. GT. N0VFL03 GO TO 106 
IF(lEX3< IK3. GT. N0VFL03 GO TO 106 
IF‘'LEX4( IM3. GT. N0VFL03 GO TO 106 


IE:V. ERRl 3 


CONTI NUE 
GO TO 108 
IEU=-S010 
RETURN TERRI 

THE LAST MM ROWS ARE COMPUTED SEPARATELY 
DO LOOP TO GO THROUGH THE MM ROWS 

DO 140 LR=1.MM 
ISR=IPT(LR+1 3 

DO LOOP TO GO THROUGH EACH ROW MM TIMES 

DO 142 IRW=1.MM 


I = IDATAn:RW. ISR3-LEAST1 

CRIOMM. FAGH) lo 
OF FO-OIi aUALFTY 


DO LOOP TO WORK DOWN A ROW, COMPUTING 
THE 0-DEGREE NEIGHBOUR ONLY 

DO 144 N=IRW, NIJMPMM, MM 
NMM=N+MM 

J= I DATA ( NMM , I SR ) -LEAST 1 


I J=INDEX< I, J) 

COUNT HORIZONTALLY ADJACENT ( 0-DEGREE ) 
NEIGHBOUR 


LEX2( IJ)=LEX2(IJ)+1 


144 I=J 
142 CONTINUE 

140 CONTINUE 


DOUBLE THE DIAGONAL TO MAKE EVERYTHING 
COME OUT RIGHT 

NOBL.= IMAX-I MIN+1 
DO 12 I=1,N0BL 

II=INDEX(I,I) 

LEX 1(11) =LEX 1 ( 1 1 ) *2 
LEX 2(11) =LEX 2 (II) -«-2 
LEX3( II )=LEX3( II )*2 
12 LEX4(II)=LEX4( II)*2 

CALL CLOSE (I DAT I) 

CALL KDPOP 
RETURN 

ERROR 

9999 CALL CLOSE ( IDATl ) 

RETURN lERRl 

END 
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F-IJ-M-C-:i 


CFLINCI 

c 

C' 

C PROGF^AM 

C TITLE 

C 

C PRni.M.:AMNER 

C UPDATE 

C: 

C: 

C 

C DnCUMEN- 

C TAT I ON 

C 

C COMPUTER 

C REQUIRED 


ij 

C 

c 


u 

c 

c 

c 

c 

c 

c 

c 

C: 

C 

C 

C 

C 

C 


PROGRAM 

language 

PURPOSE 


METHOD 


CALLING 

SEGliJENCE 

ARGUMENTS 




PARAMETERS 
AND ARRAYS 


INPUT AtMD 
OUTPUT 

RETURNS 


SUBROUTINE FUNGI 


A. SINGH OCTOBER 1972 
ROBERT M HARALICF-.: FEDRUARY 1974 

GE MONAGHAN AUGUST S. 1974 

CH IN-HUANG CHEN FEBRUARY 22. 1975 

A. SINGH 


ANY 


FORTRAN lY 


FUNCl COMPUTES THE SUM PROBABILITY FEATURE OF THE. 
IMAGE. 

FUNCl FIRST COMPiJTES THE TOTAL NUMBER OF PAIRS FOR 
EACH DIRECTION. THEN P(I,J) FOR THE K! LEX ARRAY IS 
P(I.. J) = LEXHU: IJ)./-<TOTAL NUMBER OF PAIRS FOR THE F=:; 
LEX^ ARRAY ) . I J •■= I NDEX (I . J ) . 

CALL FUNCl (LEX 1. LEX2. LEX3, LEX4. FUNC, NBUBL) 


LEXi. LEX2. LEX3 AND LEX4 ARE THE FOUR TRIANGULAR 
GRAY TONE MATRICES. 

FUNC THIS IS A TWO DIMENSIONAL ARRAY WHERE THE 

RESULTS OF SUBROUTINE FUNCl ARE STORED. 
THESE ARE STORED IN TRIANGULAR FORM LIKE 
THE LEX ARRAYS. THE SECOND SUBSCRIPT 
CORRESPONDS TO THE DIRECTION ( H:'.= l .. 2.. 3 OR 4 
IS 90.0.135 OR 45 DEGREES RESPECT I YELY ) . 
WHILE THE FIRST SUBSCRIPT. I J=INDEX ( I , J ) , 

IS THE LOCATION OF THE GRAY TONE PAIR (I.J) 
AS IN THE LEX ARRAYS. 

NBUBL SIZE OF A LEX ARRAY 

NBUBL=NOBL-K- ( NOBL+ 1 ) /2 

NOEL NUMBER OF GRAY TONES 

Ri.. R.2.. R3. R4 ARE THE. REC I PR I CAL OF THE TOTAL NUMBER 
OF GRAY TONE PAIRS FOR EACH OF THE FOUR 
DIRECTIONS. 

NONE 


WO ERROR RETURNS 
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INDEX 


SUBPROGRAMS 
REQUIRED 

CALLED BY TXTMN 

SUBROUTINE FUNCi ( LEX 1 , LEX2. LEX3> LEX4, FUNC. NBUBL. NOBL) 
DOUBLE INTEGER FUNC, LEX 1 , LEX2, LEX3, LEX4 

DIMENSION LEXK i ). LEX2( 1 ) , LEXSd), LEX4( 1 ). FUNC( 4) 

FUNi::(NBUBL, 4) 


NOW COMPUTE FUNC 

TO DETERMINE THE TOTAL NUMBER OF PAIRS IN A GIVEN DIRECTION 

Ri=0. 

R2=U. 

R3=0. 

R4=0. 

DO 5 1=1, NOBL 
DO 5 J=l, NOBL 
I J=INDEX( I, J) 

TEMP=LEX1 ( IJ) 

R1=R1+TEMP 
TEMP=LEX2( IJ) 

R2=R2+TEMP 
TEMP=LEX3< IJ) 

R3=R3+TEMP 
TEMP=LEX4( IJ) 

R4=R4+TEMP 

CONTINUE 

TO COMPUTE AVERAGE 

AVG1=0. 

AVG2=0, 

AVG3=0. 

AVG4=0. 

DO 6 1=1, NOBL 
DO 6 J=l,NOBL 
IJ=INDEX ( I, J) 

TEMP=LEX1 •; IJ) 

AVG 1=AVG 1 +TEMP-K-TEMP 
TEMP=LEX2< IJ) 

AVG2=AVG2+TEMP*TEMP 
TEMP=LEX3< IJ) 

AVG3=AVG3+TEMP’«-TEMP 
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TEMP - I EX 4 ( ;rj j 

AVG4^AVG4+TEMF-i!-rEMP 
C014T T Ml IE 


A V G I ~ A V L'> 1 y R 1 
AVG2=AVG2/R2 
Pi V G •:.! A V G G / R ■«> 
A V G 4 •— A V Li 4 / R 4 


DO 7 1 = 1; NOBL 
DO 7 J=I.;NOBL 
I J= INDEX < I.. J) 

TEMP=LEX1 •; IJ) 

FUNO ( I J , 1 ) = ( TEMP- A VO 1 ) 1 000. /R I 
TEMP=LEX2': TJ) 

FUNC I J , 2 } = < TEMP-A VG2 ) 1 000, ,/’R2 
TEMP=LEX3i; IJ) 

FUMC (I J; 3 ) = ( TEMP-AVG3 ) * 1 000, ,4R3 
TEMP=LEX4< IJ) 

FUNC ( I J; 4 ) = ( TEMP- AVG4 ) 1 000. ,/R4 
CONTINUE 


RETURN 

END 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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F-U-N-C-2 


PROGRAM 

TITLE 

PROGRAMMER 

UPDATE 


DOCUMEN- 

TATION 

COMPUTER 

REQUIRED 

PRUURi-il''l 

LANGUAGE 

PURPOSE 


METHOD 


CALLING 

SEQUENCE 

ARGUMENTS 


PARAMETERS 
AND ARRAYS 


SUBROUTINE FUNC2 


A. SINGH OCTOBER 1972 

ROBERT M HARALICK FEBRUARY 1974 

GE MONAGHAN OCTOBER 1974 

CH IN-HUANG CHEN FEBRUARY 22/ 1975 

A, SINgH 


ANY 


FORTRAN IV 


FUNC2 COMPUTES THE GRADIENT ENTROPY FEATURE OF THE 
IMAGE. 

FUNC2 FIRST COMPUTES THE TOTAL NUMBER OF PAIRS FOR 
EACH DIRECTION. THE GRADIENT ENTROPY COMPONENT IS 
ALOG ( 1 . +ABS >: I - J ) ) ALOG ( P (I , J ) ) / WHERE THE PROBAB I L I TY 
IS P<I/J) = LEXKaJ)./<TOTAL NUMBER OF PAIRS FOR THE 
K LEX ARRAY ) . I J - I NDEX ( I / J ) . 

CALL FUNC2(LEXL LEX2/ LEX3/ LEX4/ FUNC/ NBUBL) 


LEXl, LEX2, LEX3 AND LEX4 ARE THE FOUR TRIANGULAR 
GRAY TONE MATRICES. 

FUNC THIS IS A TWO DIMENSIONAL ARRAY WHERE THE 

RESULTS OF SUBROUTINE FUNC2 ARE STORED. 
THESE ARE STORED IN TRIANGULAR FORM LIKE 
THE LEX ARRAYS. THE SECOND SUBSCRIPT 
CORRESPONDS TO THE DIRECTION (K=L2l3 OR 4 
IS 90/0/135 OR 45 DEGREES RESPECTIVELY)/ 
WHILE THE FIRST SUBSCRIPT/ I J=INDEX ( I / J ) / 

IS THE LOCATION OF THE GRAY TONE PAIR <I/J) 
AS IN THE LEX ARRAYS. 

NBUBL SIZE OF A LEX ARRAY 

NELlBL=NOBL-«- C WOBL+l ) /2 

NOBL NUMBER OF GRAY TONES 

R1/R2/R3/R4 ARE THE RECIPRICAL OF THE TOTAL NUMBER 
OF GRAY TONE PAIRS FOR EACH OF THE FOUR 
DIRECTIONS. 

RL 1 .. 2/ 3 / 4 ARE THE PROBAB I L I T I ES P ( I / J ) / FOR THE 

FOUR DIRECTIONS/ FOR GRAY TONE I TO OCCUR 
NEXT TO GRAY TONE J IN A PARTICULAR 
DIRECTION. 
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INPUT Al'lD 
OUTPU'T 

RETURNS 

SI IBPROGRAMS 
REQUIRED 

CAL.LED BY 


MO ERROR RETURNS 
I MDEX 

TXTMN 


SUBROUTINE FUNC 2 ( LEX 1 . LEX 2 , LEX 3 , LEX 4 , FUNC, NEUBL. MOBL) 

DOUBLE INTEGER FUNO, I..EX 1 / LEX 2 . LEX 3 . LEX 4 
D I MENS I ON LEX 1 ( 1 ), LEX 2 ( 1 ) , LEX 3 ( 1 ) , LEX 4 ( 1 ), FUNC ( 1 , 4 ) 

FUNCRNBUBL, 4 ) 


fVl M. 

■ . 

R3^=-0. 

R 4 =^ 0 . 


NOW COMPUTE FUNC 

TO DETERMINE THE. TOTAL NUMBER OF PAIRS IN A GIVEN DTPPrTlCiM 


DO 5 I=U. ,NOBL 
DO '~i J=l,NOBL 
I J= INDEX ( I, U; 

TEMP“LEX 1 ( I J) 
Rl-^Rl+TEMP 
TEMP=LEX 2 ( IJ) 
R 2 "R 2 +rEMP 
TFMP=LEX 3 ( IJ) 
R 3 =R 3 +TEMP 
TEMP=LEX 4 ( I.-,') 
R 4 =R 4 +TEMP 
5 CONTINUE 


TO GET RT, R 2 , R 3 , R 4 TO SAVE DIVISION'S 

R.l-t. ,/Rl 
R 2 ~l /R 2 
R .1— 1 ’ R 
R 4 ~.l, ,/R 4 
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TO COMPUTE ANGULAR MOMENTUM COMPONENT 


DO 2 I=i,NOBL 
DO 2 J =1, NOBL 
I J= INDEX (I. J) 

TEMP=LEX1 ( IJ) 

RL1=TEMP#R1 
TEMF-LEX2a J) 

RL2=TEMP*R2 
TEMP=LEX3( IJ) 

RL3=TEMP*R3 
TEMP--=LEX4( IJ) 

RL4=TEMP«-R4 

FIJNC ( I J, 1 ) -ALOG ( 1. +ABS < FLOAT ( I-J ) ) ) *ALOG < 1. E-9+RL1 ) -«-200. 
FUNC< IJ. 2)=AL0G( 1. +ABS( FLOAT ( I-J) ) )->^ALOG< 1. E-9+RL2 ) . 
FIJNC ( I J> 3 ) =ALOG ( 1 . +ABS < FLOAT (I-J ) ) ) *ALOG ( 1 . E-9+RL3 ) #200. 
FUNi::< IJ. 4)=AL0G( 1. +ABS( FLOAT < I-J) ) )#ALOG< 1. E-9+RL4 ) #200. 

2 CONTINUE 

RETURN 

END 
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PtVi ir^PF,M 
riTLF 

PFVrir,Rf:,M|v|£;p 

UPDATE 


DOCUMEN- 

TATION 

COMPUTEFv 
REQU 1 RED 

PP!OOF'‘hI''I 
L AiNUUhOE. 

PURPOSE 

METHOD 


CALI... IMG 
SEQUENCE 

ARGUMENTS 


E'ARAME TERsi 
AMD ARRAYS 


INPUT AMD 


SUF-'POUTIME FLINCH 


A. SINGH OCTOBER 1972 
ROBERT M HARALICK FEBRUARY 1974 
GE MONAGHAN OCTOBER 1974 

CHTN-HUANG CHEN FEBRUARY 22. 1975 


A. SINGH 


AMY 

FORTRAN IV 


FUNCS COMPUTES THE ENTROPY FEATURE OF THE IMAGE. 

FUNC3 FIRST COMPUTES THE TOTAL NUMBER OF PAIRS FOR 
EACH DIRECTION. THE ENTROPY COM’^’ONEMT IS THEN 
- P ( I . J ) -frALOG >: P ( I WHERE THE PROBAE I L I TY P ( I . J ) 

IS pa..j) = LEXK.';rj)/rroTAL number of pairs for 
THE K L_EX ARRAY). IJ - INDEX (I..U) 

CALL FUNC3(LEX1. LEX2. LEX3. LEX4. FUNC. NBUBL) 


LEXl. LEX2. LEX3 AND LEX4 ARE THE FOUR TRIANGULAR 
GRAY TONE MATRICES. 

FUNC THIS IS A TWO DIMENSIONAL ARRAY WHERE THE 

RESULTS OF SUBROUTINE FUNC3 ARE STORED. 
THESE ARE STORED IN TRIANGULAR FORM LIKE 
THE LEX ARRAYS. THE SECOND SUBSCRIPT 
CORRESPONDS TO THE DIRECTION (K=1.2. 3 OR 4 
IS •-'0,0.135 OR 45 DEGREES RESPECTIVELY). 
WHILE THE FIRST SUBSCRIPT. I J=INDEX (I. U) . 

IS THE LOCATION OF THE GRAY TONE PAIR (I.J) 
AS IN THE LEX ARRAYS, 

NBUBL SIZE OF A LEX ARRAY 

NB1.JBL=N0BL*- ( NOBL+ 1 ) /2 

NOBL NUMBER OF GRAY TONES 

Pi.R2.R3. R4 ARE THE RECIPRICAL THE TOTAL NUMBER 
OF GRAY TONE PAIRS FOR EACH OF THE FOUR 
DIRECTIONS. 

RL1.2. 3. 4 ARE THE PROBABILITIES Pa. J) . FOR THE 

FOUR DIRECTIONS. FOR GRAY TONE I TO OCCUR 
NEXT TO GRAY TONE J IN A PARTICULAR 
DIRECTION, 


ORIGINAL PAGB IS 
OP P(X)R QUALilY 


175 


OUTPUT 


RETURNS NO ERROR RETURNS 

SUBPROGRAMS I NDEX 
REQUIRED 

CALLED BY TXTMN 

SUBROUTINE FUNCSKLEXl. LEX2, LEX3, LEX4, FUNC, NBUBL. NOBL) 

DOUBLE INTEGER FUNC, LEX 1 . LEX2. LEX3, LEX4 
D I MENSION LEX 1< 1 ) , LEX2 ( 1 ) , LEX3 ( 1 ) , LEX4 ( 1 ) . FUNC (1.4) 

FUNC(NBUBL, 4) 


NOW COMPUTE FUNC 

TO DETERMINE THE TOTAL NUMBER OF PAIRS IN A GIVEN DIRECTION 


R1=0. 

R2=0. 

R3=0. 

R4=0. 

DO 5 1=1. NOEL 
DO 5 J=l,NOBL 
IJ=INDEX( I. J) 
TEMP=LEX1 < IJ) 
R1=R1+TEMP 
TEMP=LEX2( IJ) 
R2=R2+TEMP 
TEMP=LEX3( IJ) 
R3=R3+TEMP 
TEMP=LEX4( IJ) 
R4=R4+TEMP 
CONTINUE 


TO GET R1.R2. R3. R4 TO SAVE DIVISIONS 


Rl = l. /R1 
R2=l. /R2 
R3=l. /R3 
R4=l. /R4 

TO COMPUTE ENTROPY COMPONENTS 

DO 2 1=1, NOBL 
DO 2 J=I.NOBL 
I J=INDEX( I. J) 


176 


TEhP=^L.E X J. ( I J) 
Rl.i TEMP-K-Rl 
I Pl’ih‘=LLX wi ( I. J; 
Rl„.? = TEMPttR2 
TEMP^^=IEX3( i: J) 
RL3 =- TEMP^^rR3 
TFMP-=L.EX4': IJ) 
RL4 - TEMP-K-R4 


IF(RLt LT, 0. 00000 IJ GO TO OH 
FIJMC ■: I .J. 1 ) == ■: -RL :! -K-ALOG < RL ;l. ) ) #200. 
IF ( RL2. LT. 0. 000001 ) GO TO 32 
FLING ( I J. 2 ) - ( -RL2#AL0G ( RL2) ) #200. 
IF^RL.G. LT. 0. 000001) GO TO 33 
FLINlH I-J/ O) = ( “RlLO#ALUO ( RLO ) ) #200. 
IF’(RL4. LT. 0. 000001) GO 7”0 2 
FIJMC < I J , 4 ) "• ( -RL4# ALOG ( RL4 ) ) #200. 
CONTINUE 

RETURN 

END 


OS' POOR 
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CFU1MC4 

C 

c 

C PROGRAM 

C TITLE 

C 

C PROGRAMMEBl 

C UPDATE 



DOCUMEN- 

TATION 

COMPUTER 

REQUIRED 

PROGRAM 

LANGUAGE 

F'URlF’UGE 


METHOD 


CALLING 

SEQUENCE 


F-U-N-C-4 


SUBROUTINE FUNC4 


ROBERT M HARALICK MAY 1973 
ROBERT M HARALICK FEBRUARY 1974 
GE MONAGHAN OCTOBER 1974 

CH IN-HUANG CHEN FEBRUARY 22, 1975 

ROBERT M HARALICK 


ANY 


FORTRAN IV 


FUNC4 COMPUTES THE GRADIENT FEATURE OF THE 
IMAGE. 

F1JNC4 FIRST COMPUTES THE TOTAL NUMBER OF PAIRS FOR 
EACH DIRECTION. THE GRADIENT COMPONENT IS 
ABS < I - J ) ,/P C I , J ) WHERE THE PROB AB I L I TY P ( I , J ) I S 
P(I,J) = LEXKaj)./(TOTAL NUMBER OF PAIRS FOR THE 
K LEX ARRAY ) . I J = I NDEX ( I , J ) . 

CALL FUNC4(LEX1, LEX2, LEX3, LEX4, FUNC, NBUBL ) 


C; ARGUMENTS LEX 1 , LEX2, LEX3 AND LEX4 ARE THE FOUR TRIANGULAR 

C GRAY TONE MATRICES. 

C: FUNC THIS IS A TWO DIMENSIONAL ARRAY WHERE THE 

C RESULTS OF SUBROUTINE FUNC4 ARE STORED. 

C THESE ARE STORED IN TRIANGULAR FORM LIKE 

C THE LEX ARRAYS. THE SECOND SUBSCRIPT 

C CORRESPONDS TO THE DIRECTION (K=l>2,3 OR 4 

C IS 9Ch 0,135 OR 45 DEGREES RESPECTIVELY), 

C WHILE THE FIRST SUBSCRIPT, I J=INDEX ( I , J ) , 

C IS THE LOCATION OF THE GRAY TONE PAIR (I,J) 

C AS IN THE LEX ARRAYS. 

C NBUBL SIZE OF A LEX ARRAY 

C NBUBL=N0BL*(N0BL+1 )./'2 

C ■ 

C PARAMETERS NOBL. NUMBER OF GRAY TONES 

C AND ARRAYS Rl, R2, R3, R4 ARE THE RECIPRICAL OF THE TOTAL NUMBER 

C OF GRAY TONE PAIRS FOR EACH OF THE FOUR 

C DIRECTIONS 

C RL1,2, 3, 4 ARE THE PROBABILITIES PM,J), FOR THE 

C FOUR DIRECTIONS, FOR GRAY TONE I TO OCCUR 

C NEXT TO GRAY TONE J IN A PARTICULAR 

C DIRECTION. 
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IMPlJl AND 
OUT PI. IT 


NONE 


RETURNS NO ERROR RETURNS 

SUBPROGRAMS IlvlDEX 
REQUIRED 

CALLED GY TXTMN 

SUBROUTINE FUNC4<LEX.l, LEX2, LEX3, LEX4. FLINCY NBUBL,. NOBL) 
DQUBL.E INTEGER FUMC, LEXl/ LEX2, LEX3, LEX4 

D 1 MENS I ON L. E X 1 (1 ) . LEX 2 ( 1 ) , LEX 3 U ) . LEX 4 (1 ) . FUNC 1 . 4 ) 

FUMC < NBUBL. 4) 

AF= 1 . /FLOAT •: NBUBL ) 


MOW COMPUTE FUNC 

TO DETERMINE THE TOTAL NUMBER OF PAIRS IN A GIVEN DIRECTION 

R1=0. 

R3"0. 

R4=---0^ 

DO 5 1=1. NOBL 
DO 5 J=l, NOBL 
I .J=INDEX ( I. J) 

TEMP=LEXiaj) 

R1=R1+TEMP 
TEMP=LEX2< IJ) 

R2=R2+TEMP 

TEMP=LEX3(IJ) 

R3=R3+TEMP 
TEMP=LEX4< I.J> 

R4=R4+TEMP 

CONTINUE 

TO GET R1.R2. R3.R4 TO SAVE DIVISIONS 

Rl = l, /R1 
R2=l. /R2 
R3=l /R3 
R4=L /R4 

TO COMPUTE ANGULAR MOMENTUM COMPONENT 

DO 2 1=1. NOBL 
DO 2 .J =1. MOBL 
IJ=INDEXa, ...n 


179 


TEMP=LEXU IJ) 

TEMP=LEX2( IJ) 

FUNC ■: IJ, 1 ) = ( ABS ( FLOAT ( I -J ) ) / ( AF+TEMP*R 1 ) ) ■t4-200. 
TEMP=LEX3( IJ) 

FUNC < I J, 2 ) = ■: ABS ( FLOAT ( I - J ) ) / ( AF+TEMP*R2 ) ) -t^200. 
TEMP=LEX4( IJ) 

FUNC ( I J, 3 ) •■= C ABS ( FLOAT ( I -.J ) ) / ( AF+TEMP-X-R3 ) ) -t«-200. 

FUNC < I J, 4 ) - ( ABS ( FLOAT ( I - J ) ) / ( AF+TEMP^^-R4 ) ) -K200. 
CONTINUE 

RETURN 

END 
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CFUMCS 


F-IJ-N-C-S 


FLEX XJ QUANT I Z I NG FUNCT I ON 


PROGRAM TITLE:' FUNC5 

VERSION; A 

DATE; NOVEMBER 23, 19'73 

AUTHOR: ROBERT M HARALICK 

UPDATE CH IN-HUANG CHEN 2/22/75 

DOCUMENTED BY: ROBERT M HARALICK 

IMPLEMENTED ON; PDP 15 

LANGUAGE: FORTRAN 

PURPOSE: 

TH I S SUBROUT I ME PREPARES NORMAL I Z ED LE'X 
ARRAYS WHICH HAVE BEEN EQUAL PROBABILITY QUANTIZED ACCORDING 
TO THEIR DIAGONAL ELEMENTS AND PUTS THE RESULTS IN FUNC ARRAY. 


ENTRY POINT; 


FLINCS'ILEXI, LEX2, LEX3, LEX4, FUNC, NBUBL) 


ARGUMENT LIST: 


LEX 1 
LEX 2 
LEX 3 
LEX 4 
FUNC 


NBUBL 


IS VEIRTICAL CO-OCCURENCE MATRIX 
I S HOR I Z ONT AL CO-OCCURENCE MATR I X 
IS 135 DEGREE CO-OCCURENCE MATRIX 
IS 45 DEGREE CO-OCCURENCE MATRIX 
IS THE NORMALIZED AND QUANTIZED 
CO-OCCURENCE MATRICES 
FUNC- NOBL, 4) 

IS THE SIZE OF THE LEX ARRAYS 


SUBROUTINE FUNCS ( LEX 1 , LEX2, LEX3, LEX4, FUNC, NBUBL, NOBL) 
DOUBLE INTEGER FUNC, LEX 1, LEX2, LEX3, LEX4, F 

D I liEMS I ON LEX 1 >; 1 ) , LEX2 < 1 '1 , LEX 3 ( 1 , LEX 4 ( 1 ) , FUNC ( 1 , 4 ) , F C 1 ) 
DATA INTVD /S/ 

CALL ADJl (F, FUNC< 1, 1 ) ) 


CALL LEXEQPa.EX4, NOBL, INTVD, F) 
DO 12 1-1, NBUBL 
FUNC< I, 4)=-FC I) 


CALL LEXEQP<LEX3, NOEU, INTVD, F) 
DO 13 1=1, NBUBL 
FUNCn,3)=F(I) 


CALL LEXEQP< LEX2, NOBL, INTVD, F) 
DO 14 1 = 1, NBUBL 
FUNC’ I, 2)=F(X) 


ORIGINAL PAGE IS 
OF POOR QUALnY 


CALL LEXEQP ( LEX 1 , NOBL, I NT VD, F ) 
RETURN 

DO 15 1=1, NBUBL 
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FIJNC< I, 1 )=-F< I ) 
END 
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CLEXEQP 

C 

c 

c 


c 

c 

c 

c 

c 

c 


L ~ E — X ~ E “ U “ P 

EQUAL PROBABILITY QUANTIZE THE DIAGONAL OF THE LEX ARRAY 

LEXEQP 


PROGRAM TITLE: 
VERSION: 

DATE: 

AUTHOR: 

UPDATE 

DOCUMENTED BY: 

LANGUAGE: 

PURPOSE: 


A 

NOVEMBER 2:5, 197:3 
ROBERT M HARALICK 
CH IN-HUANG CHEN 2/22/7!: 
ROBERT M HARALICK 
FORTRAN IV 


THIS SUBROUTINE EQUAL PROBABILITY QUANTIZES 
THE LEX ARRAY ON THE BASIS OF THE DIAGONAL ELEMENTS. 


ENTRY POINT: LEXEQP (LEX, NOBL, I NTVD, FUNC) 

ARGUMENT LIST: 


LEX 

NOBL 

INTVD 


FUNC 


IS THE LEX ARRAY 

IS THE NUMBER OF BRIGHTNESS LEVELS 
IS THE NUMBER OF DESIRED QUANTIZED 
LEVELS 

IS THE NORMALIZED AND QUANTIZED LEX 
ARRAY. 


SUBROUTINE LEXEQP ( LEX > NOBL, INTVD, FUNC) 
DOUBLE INTEGER FUNG, LEX 


C 

_ C 
( 


DIMENSION LEX ( 1 ) , FUNC< 1 ) 

COMMON / 1 0/NS I ZE, F ( 1 6 ) , FLQ < 1 6 ) , MEX (1 :3S ) , I T ( 1 76 ) 


PUT CUMULATIVE DISTRIBUTION OF 
DIAGONAL ELEMENTS OF LEX ARRAY INTO F. 


I F ( I NT VD. GT, 16) I NT VD= 1 6 
NSIZE=1024 

NBUEL=NOBL«- ( NOBL+ 1 ) /2 
S=0 ' ' 

DO 1 1=1, NOBL 
II=INDEX< I, I ) 
TEMP=LEX(II) 

1 S=S+TEMP 

s=i . /s 

TEMP=LEX ( 1 ) 

F(1)=TEMP*S 
DO 2 I=2,N0BL 
J=INDEX( I, I) 

TEMP=LEX(J) 

2 F ( I ) =F ( I - 1 ) +TEMP-K-S 
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I c 
i c 
I C 
^ c 


12 



c 


13 




10 

7 

4 



RCiBL.=FLOAT ( NOBL) 

CALL EQPQNKNOBL, INTVD, F, FLQ, ROBL. 0, 01 ) 


CONSTRUCT THE QUANTIZED LEX MATRIX 

Jl = l 

DO 4 J=l, INTVD 
IF( J. EQ. 1 ) GO TO 12 
J1=FLQ( J-1 ) + l. 

CONTINUE 

J2=FLQ(J) 

Kl = l 

DO 7 K=1, J 

IF(K. EQ. 1 ) GO TO 13 

K1=FLQ<K-1 )+l. 

CONTINUE 

K2=FLQ(K) 

MM=INDEX<: J, K) 

MEX(MM)=0 

DO 10 JJ=J1, J2 
DO 10 KK==K1.K2 
LL=INDEX( JJ. KK) 

MEX ( MM ) =MEX < MM ) +LEX ( LL ) 

CONTINUE 

CONTINUE 

CONTINUE 

DEFINE THE QUANTIZING FUNCTION 

J=1 

DO 3 I=1.N0BL 

IF(FLOAT< I). LE. FLQ< J) ) GO TO 5 

GREY TONE I BELONGS TO THE NEXT QUANTIZING 
INTERVAL, 


J=J+1 

GREY TONE I EsELONGS TO THE JTH QUANTIZING 
INTERVAL. 

IT< I ) -• J 
CONTINUE 

TRANSFER IT TO FUNC. 

DO 11 1 = 1, NOBL 

II=IT(I ) ■ . 

DO 11 J=l, NOBL 
JJ=IT'; J) 

N=INDEXM, J) 

MM=INDEX< ri, JJ) 

TEMP=MEX(MM) 
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CONTINUE 

FUNC ( IM ) =TEMP*S« 1 000. 

RETURN 

END 
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P-L-X-l-T 


CPLXIT 
C 

c 

C PROGRAM SUBROUTINE PLXIT 

C TITLE 


C 

C 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAMMER A. SINGH NOVEMBER ?<72 

MODIFIED 5/14/73 ROBERT M HARALICK 

7/10/73 

2/2/74 

8/10/74 GE MONAGHAN 
10/10/74 RM HARALICK 
2/22/75 CH IN-HUANG CHEN 


DOCUMEN- A. SINGH 

TATION 

COMPUTER PDP 15 

REQUIRED 

PROGRAM FORTRAN IV 

LANGUAGE 


PURPOSE PLXIT COMPUTES THE JDATA IMAGE 


METHOD PLXIT COMPUTES THE JDATA IMAGE UTILISING THE RESULTS 

OF FPLXIT AND FUNC. LET G(I,J) BE THE GRAY LEVEL OF 
THE JTH RESOLUTION CELL IN THE ITH LINE OF THE 
CONSIDERED IMAGE < I DATA ); AND LET V<I,J) BE THE JTH 
RESOLUTION CELL IN THE ITH LINE OF THE JDATA IMAGE. 
THEN — 

V < I . J ) = FUNC < G ( I , J+L ) , G < I -L, J+L ) , G ( I ~L. J ) , G ( I ~L, J-L ), 

GC I, J-L) , G< I+L. J-L) , G< I+L. J) , G< I+L, J+L) ) , 



C 

C 

c 

c 

c 


c 


i; 




WHERE FUNC IS A FUNCTION (SUCH AS FUNCl, FUNC2 OR 
FUNC3) PROVIDED BY THE USER. 

L = 1, 2. 3. > IS THE 

SEPARATION BETWEEN CELLS. L=l, MEANS NEAREST 
NEIGHBOURS L=2. MEANS NEXT TO NEAREST NEIGHBOUR ETC 
PLXIT WORKS FOR ALL POSITIVE L. 

ENTRY POINT PLXIT < IDATI. IDATJ, IDATA, JDATA, IDENT; FUNC, IPT 

NBUBL, MMl, NXMIN, NXMAX, JDENT, lEV, lERRl ) 

ARGUMENTS IDATI DAT SLOT WHERE ORIGINAL IMAGE RESIDES 

IDATJ DAT SLOT FOR JDATA IMAGE 

IDATA SCRATCH ARRAY WHERE THE ORIGINAL IMAGE IS 

READ IN. 

JDATA INTEGER ARRAY WHERE THE JDATA IMAGE IS 

GENERATED AND STORED BEFORE BEING WRITTEN 
ONTO THE TAPE (03). 

IDENT IDENTIFICATION ARRAY OF IDATA 

JDENT IDENTIFICATION ARRAY OF JDATA 

FUNC A TWO DIMENSION ARRAY CONTAINING THE 


RE’-.l ll TS OF THE EXTFRIvlA! FI HMCTinN PROGRAM. 

■ THE SECOMO- INDEX DETERMINES THE DIRECTION; 
WHILE THE FIRST ONE CGRRESPCiMDS TO THE 
ELEMENT IN THE ASSOCIATED LEX ARRAY. 

IPT ARRAY WHICH CONTAINS THE POINTERS FOR 

THE IDATA AMD THE JDATA ARRAYS 
NBLiBL SIZE OF A LEX ARRAY 

NBUBL=MOBL#(NOr;!L+l WHERE NOEL IS THE 

NUMBER OF GRAY TONES. 

MMl SPATIAL DISTANCE +1 

NX MIN MINIMUM JDATA VALUE 

NX MAX MAXIMUM JDATA VALUE 

lEV INTEGER EVENT VARIABLE 

IEV=-5011 IF NUMPPL OR NUMLIN IS LESS THAN 
TWICE MM. 


PARAMETERS 


AND ARRAYS 


NUMLIN NUMBER OF LINES IN THE IMAGE 

NUMPPL NUMBER OF POINTS PER LINE IN THE INPUT IMAGE 

IMAX LARGEST GRAY TONE ON INPUT FILE 

IMIN LEAST GRAY TONE ON ILJF'UT FILE 

LEAST 1 =IMIN-1. LEAST! IS USED FOR NORMALIZING 
THE GRAY TONES. 

NOBL NUMBER OF GRAY TONES 

NOBL=IMAX-IMIN+l 


RETURNS 


NO ERROR RETURNS 


SUBPROGRAMS INDEX 
REQUIRED 

CALLED BY TXTMN 

COMMENTS PLX IT WORKS FOR ALL SPATIAL DISTANCES, IT DOES 

THIS E!Y HAVING MM + 1 LINES OF IDATA IN THE CORE; 
WHERE MM IS THE SPATIAL DISTANCE PARAMETER. 

SUBROUTINE PLXITCIDATI; IDATJ.- IDATA; JDATA; IDENT; FUNC; IPT; NBUBL; MMl. 
i IMGNO; NXMIN; NXMAX; JDENT; lARRAY; lEV; lERRl ) 


DOUBLE INTEGER FUNC 

DIMENSION IDATA( 1; 1 ) ; JDATA*: 1; 1 ); FUNC( 1; 1 ) ; IPT( 1 ) ; IDENT<20) 
DIMENSION JDENT (20); I ARRAY ( 1 .. 1 , 1 ) 

IDATA (NUMPPL; MMl ) ; JDATA (NUMPPL; MM! ) ; FUNC( NBUBL; 4) ; IPT( MMl ) 


STACK SUBROUTINE NAME IN ERROR STACK 
CALL KDPUSH( •■pLXXT-'; •- ■- 3 


SET parameters 




IS 
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100 


111 

110 


NUMPPL=IDENT(6) 
NUMLIM=IDENT<7) 
IMIN=IDEMT< 15) 
IMAX=IDENT< 16) 

LEAST1 = IMIN--1 
M0niL= 1 MAX “LEAST 1 
NBUBL=N0BL* ( NOBL+1 ) /2 


MX MI N= 13 1000 
NXMAX=-131000 

CHECK IF SIZE OF IMAGE IS TOO SMALL, 
RELATIVE TO THE SPATIAL DISTANCE 
PARAMETER. 

MM=MM1-1 

MM2=MM*-2 

IF(NUMPPL. LT. MM2. OR. NUMLIN, LT. MM2) G0T0 99SS 


ZERO OUT THE JDATA ARRAY 

N1JMPMM=NUMPPL-MM 
NUMLMM=NUMLIN-MM 
DO 100 1=1, MUMPPL 
DO 100 J=1,MM1 
JDATA<I,J)=0 


READ IN THE FIRST MMl LINES OF THE IMAGE 
AND SET UP POINTERS 

DO 110 IY=1, MMl 
IPT(IY)=IY 

CALL RREADdDATI, lARRAY, IMGNO, lY, 1, IDENT, IEV,ERR1) 

DO til LY=1,NUMPPL 
IDATA(LY, IY)=IARRAY(t, 1,LY) 

CONTINUE 


SETTING UP POINTERS FOR THE FIRST AND 
LAST ROWS OF THE IMAGE ARRAYS 

IST=IPT< 1 ) 

LST=IPT<MM1) 

GO THROUGH ALL BUT MM ROWS OF IMAGE 

NEXT=MM1+1 

DO 105 LCNT=1, NUMLIN 

GO THROUGH EACH ROW MM TIMES. THE FIRST 
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SET GF MM COLUMNS ARE HANDL.ED SEPARATELY 


DO 120 ]:RW=U.MM 
IRM=:IRW+MM 

SET I, L, J AMD K EQUAL TO THE 
(NORMALISED) VALUES OF GRAY TONES OF 
RESOLUTION CELLS IN POSITIONS A, B. C 
AND D AS IN THE DIAGRAM — 

A C 
B D 

WHERE A INITIALLY IS THE UPPER LEFT 
CORNER CELL. THE CELLS ARE A DISTANCE 
MM APART. 

I = IDATA(IRW. IST)-LEAST1 
L= I DATA ■: I RW.. LST ) -LEAST 1 
K=IDATA( IRM.LST) -LEAST! 

U=IDATA(IRM, IST)-LEAST1 


PUT THE TWO DIMENSIONAL INFORMATION 
INTO ONE DIMENSIONAL FORM. THE FUNCTION 
NEEDED TO CONVERT A DOUBLE SUBSCRIPTED 
ARRAY,. I MM ( X .■ Y ) . I NTO A S I NGLE 
SUBSCRIPTED ARRAY. IMM(Z). IS OF THE 
FORM G(X) + F(Y). WHERE G(X ) - < X-1 ) ■>«-X./2 
AND F(Y) = Y. THEREFORE 
Z = (X-1. )-fi-X./2. 

THIS IS DONE IN THE PROGRAM BY THE 
EXTERNAL FUNCTION INDEX (X.Y). 


SINCE THE ORDER OF OCCURRENCE OF THE 
GRAY TONES BELONGING TO A RESOLUTION 
CELL PAIR IS IMMATERIAL. THE ARRAYS ARE 
SYMMETRIC. WE LET THE LARGER OF THE TWO 
HAVE THE FIRST SUBSCRIPT^ I. E. , THE ARRAY 
IS STORED IN LOWER TRIANGULAR FORM. THE 
ORDER OF THE SUBSCRIPTING IS AS FOLLOWS — 
1MM( IV 1 ) = IMM( 1 ) , 

IMM(2. 1) = IMM<2) , 

IMM<2. 2) = IMM(3) . 

IMM<3. 1) IMM<4). 


I MM •: NOBL. NOBL ) - I MM ( NBUBL ) . 



THE SCANNING PROCEDURE. THAT IS THE 
METHOD BY WHICH THE PAIRWISE COMPARISONS 
ARE MADE. IS DESCRIBED BELOW FOR THE 
GENERAL CASE. 

CGNSIDER A RESOLUTION CELL WITH SPATIAL 
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IL=IWDEX( I. L) 


COORDINATES <M.N), AND CALL THIS CELL I. 
THE SCANNING OPERATION BEGINS IN THE 
UPPER LEFT HAND CORNER OF THE IMAGE AND 
IT THEN PROCEEDS BY COMPARING THE GRAY 
TONE OF WITH AT MOST FOUR GRAY TONES 

OF ITS NEIGHBOURING RESOLUTION CELLS. 

THAT NEVER NEEDS TO CONSIDER MORE 

THAN FOUR NEIGHBOURS CAN BE SEEN FROM 
THE DIAGRAM OF THE SEARCH PATTERN SHOWN 
BELOW — 

I J 
M L K 

ON A GIVEN ITERATION, WILL LOOK FIRST 

AT ITS VERTICAL NEIGHBOUR (?/.L8<), NEXT 
AT ITS HORIZONTAL NEIGHBOUR THIRD 

AT ITS LOWER RIGHT NEIGHBOUR AND 

FOURTH AT ITS LOWER LEFT DIAGONAL 
NEIGHBOUR (?<MS<). THEN MOVES INTO 

THE POSITION OF THE LEFT-MOST RESOLUTION 
CELL OF THE PREVIOUSLY SCANNED SECOND 
ROW (THE POSITION OCCUPIED BY -?<M2^), 

THE OPERATION IS REPEATED UNTIL ALL 
NEIGHBOURING PAIRS OF RESOLUTION CELLS 
HAVE BEEN EXAMINED. THE PROCEDURE IS 
FURTHER REPEATED FOR CELLS SKIPPED OVER 
IF THE SPATIAL DISTANCE IS GREATER THAN 
ONE/ TILL ALL CELLS HAVE BEEN EXHAUSTED. 


ADD FUNC(IL/1) TO CENTER CELL AND 
90-DEGREE NEIGHBOUR 

JDATA( IRW/ 1ST) = JDATA( IRW, 1ST) + FUNC( IL, 1 ) 

JDATAURW, LST) = JDATA< IRW, LST ) +FUNC(IL/1) 

IJ=INDEX( 1/ j) 

ADD FUNC(IJ,2) TO CENTER CELL AND 
0-DEGREE NEIGHBOUR 


JD AT A ( I RW, I ST ) = JD AT A ( I RW, I ST ) + FUNC ( I J , 2 ) 
JDATAdRM, 1ST) = JDATAdRM, 1ST) + FUNCdJ/ 2) 


IK-INDEXd/ K) 

ADD FUNCdK/S) TO CENTER CELL AND 
135-DEGREE NEIGHBOUR 

JDATACIRW/ 1ST) = JDATAdRW, 1ST) + FUNCdK, 3> . 

JDATAC IRM/ LST) = JDATA( IRM, LST) + FUNCC IK, 3) 
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MI=1RW 



MOW ITERATE DOWN THE R!OW 

DO 130 W“IR1'-T NUMRMM,. MM 

NMM=N+Mr'l 

WNM-W-MM 

NI=M 

I 

M“L 
L ~k: 

J= I DATA ( MMM; X ST ) -LEAST 1 
K= I DATA •; MMM, L3T) -LEAST 1 

lL = :tMDEX(I..L) 

ADD FUNCC IL> 1 ) TO CENTER CELL AND 
90-DEGREE NEIGHBOUR 

JDATA( N, 1ST) = JDATA( N, 1ST) + FUNC C XL. 1 ) 

UDATA-; N.LST) - JDATA( M, LST) + FUNC<IL;1) 


IJ-TNDEX( X. J) 


ADD FUNCn,Jr2) TO CENTER CELL AND 
0-DEGREE NEIGHBOUR 


JDATA*: N. 1ST) = JDATA< N. 1ST) + FUNC< X J. 2) 
jDATA(NMM/ XST) - JDATA (NMM. 1ST) + FUNCdJ, 2) 


IK-=XNDEX ( XdO 

ADD FUNC(IK,S) TO CENTER CELL AND 
135-DEGREE NEIGHBOUR 

JDATA( N, 1ST) ~ JDATA ( N, 1ST) + FUNGdKS 3) 

JDATAXNMM.. LST) = JDATA ( NMM, LST ) + FUNCvIKSS) 

IM=IWDEXd^M) 

ADD FUNCC I ML4) TO CENTER CELL AND 
45-DEGREE NEIGHBOUR 

JDATA ( N.' 1ST) = JDATA( N/IST) + FL1NC< IM, 4) 

JDATA ( NNM..LST) = JDATA( NNM. LST) + FUNCdM. 4) 

130 CONTINUE 

COMPUTE THE LAS:T SET OF MM COLUMNS 
SEPARATELY 

. nin^ni+mm 

M~L- 
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L=K 

IL=IMDE>: (I, L) 

ADD FUNC<IL, i) TO CENTER CELL AND 
90-DEGREE NE I GHBOUR 

JDATA(NIM, 1ST) = JDATA(NIM, 1ST) +FUNC(IL>1) 

JDATA<NIM. LST) = JDATA(NIN. LST) + FUNC< IL, 1 ) 

IM=INDEX( I, M) 

ADD FUNCdM, 4) TO CENTER CELL AND 
45-DEGREE NEIGHBOUR 


JDATACNIM; 1ST) = JDATA< NIM, 1ST ) + FUNCdM, 4) 
JD AT A ( N I .. LST ) = JD AT A < N I ; LST ) + FUNC d M, 4 ) 

K C 

120 CONTINUE 


. C TO WRITE OUT THE COMPLETED LINE OF THE 

C JDATA IMAGE 

: C 

- DO 699 J-l.MM 

IXM=NUMPPL-J+1 

JDATA ( Jd ST ) = ( JDATAX Jd ST ) *8 ) /5 
.. JDATAdXMdST) = ( JDATA< IXMdST)-«-8)/5 

699 CONTINUE 

IF<LCNT, NE, 1) GO TO 695 

U. C 

i DO 694 J=1,NUMPPL 

I lARRAYd J) = ( JDATA< JdST)-«-5)./3 

j. 694 CONTINUE 

I GO TO 798 

i C „ 

j. 695 CONTINUE 

I DO 797 J=l. NUMPPL 

; IARRAY< Id.; J)= JDATA (J> 1ST) 

L 797 CONTINUE 

j C . , : , , 

I 798 CONTINUE 

LINE=LCNT-MM1 

i : C 

L CALL RWRITEdDATJ; lARRAYd, LINEd, JDENTdEV; ERRl ) 

I DO 700 IXM=1; NUMPMM 

j IF( JDATA< IXM, 1ST). LT. NXMIN) NXMIN==JDATAdXM> 1ST) 

IF(JDATA<IXM, 1ST). GT, NXMAX) NXMAX=JDATA<IXM, 1ST) 

I 700 CONTINUE 

I. C . 

L C SHIFT THE POINTERS FOR THE TWO ARRAYS. 

! C THIS IS DONE BY A CYCLIC ROTATION. 

I C THE POINTER ARRAY IPT IS SUCH THAT AT ANY 

C TIME THE ITH LOCATION OF IPT CONTAINS 

! C THE POINTER TO THEITH POSITION OF THE 



LINE IN I DATA OR JDATA ARRAY. FOR 
EXAMPLE. IF IPT(2)=4 THEN THE FOURTH LINE 
OV THE PHYSICAL JDATA ARRAY IS ACTUALLY 
THE SECOND LINE. AT THAT MOMENT, 


IFCLCNT EQ. NUMLIN) GO TO 105 

ROTATE IN A CYCLIC MANNER 

ITEMP^=IPT( 1 ) 

DO 135 IB=i.MM 
135 TPT(IG)=IPT( IB+i ) 

I PT( MM 1)^1 TEMP 

SET UP THE POINTERS TO THE FIRST AND 
LAST RONS OF THE TWO li'IAGE ARRAYS 

IST=IPT(i) 

LST=IPT(MM1) 

READ IN A NEW LINE INTO THE. I DATA ARRAY 

CALL rvREAD''. IDATI.. I ARRAY.. IMGNO. LCNT,. 1. IDENT. lEU, ERR!) 

DO 112 LY=1.NUMPPL 
1 12 I DATA ( LY. LST ) I ARRAY (1 . 1 , LY ) 

ZERO OUT THE LAST LINE OF THE JDATA ARRAY 

DO 145 .„ij=l. MUMPPL 
145 JDATAX JJ,. LST)=U 

105 CONTINUE 

THE LAST MM ROWS ARE COMPUTED SEPARATELY 
DO LOOP TO GO THROUGH THE MM ROWS 

ILIME=LINE 
DO 140 LR=1. MM 
ISR=IPT<LR+1 T- 

DO LOOP TO GO THROUGH EACH ROW MM TIMES 

DO 142 IRW=i.. MM 
Ir^IDATA': IRWCISR)-LEAST1 

DO LOOP TO WORK DOWN A ROW. COMPUTING 
THE 0-DEGREE NEIGHBOUR ONLY 

DO 144 N=IRW. NUMPMM. MM 
WMM=M+MM 

J=IDATAHMMM. ISR):-LEAST1 


I J= INDEX ! I. J) 
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- c 


JDATA< ^4, I SR ) 
JDATA<NMM.. I SR) 

144 I-J 
142 CONTINUE 


ADD FUNCdJ, 2) TO CENTER CELL AND 
0-DEGREE NEIGHBOUR 

JDATA( N, ISR) + FLINCaJ. 2) 

JDATA ( NMM, I SR ) + FUNC ( I J , 2 ) 


WRITE OUT THE COMPLETED JDATA LINE 


ij 

C 

C 

G 


C 


DO 698 J=1,MM 
IXM=NUMPPL-J+1 

JDATA ( J, I SR ) = ( JDATA ( J, I SR ) *8 ) /5 
JDATA( I'XM; rSR) = ( JDATAC IXM/ISR)-»-8)/5 
698 CONTINUE 


IF<LR. NE. MM) GO TO 670 

DO 696 J=1,NUMPPL 
I ARRAY < 1 , 1 , J ) - ( JDATA ( J , I SR ) #5 ) /: 
696 CONTINUE 
GO TO 896 

670 CONTINUE 


DO 897 J=1.NUMPPL 
I ARRAY < 1,1. J)=JDATA‘: J, I SR) 
397 CONTINUE 


896 CONTINUE 


- C 


701 

140 CONTINUE 


LINE-ILINE+LR 

CALL RWRITE( IDATJ, I ARRAY, 1 , LINE, 1 , JDENT, lEV.ERRl) 
DO 701 IXM=1,NUMPMM 

IF( JDATA( IXM, ISR), LT. NXMIN) NXMIN=JDATA< IXM, ISR) 
IF( JDATA-: IXM, ISR). GT. MXMAX) NXMAX=JDATA < I XN, ISR ) 
CONTINUE 


CALL ULiJSE( IDATJ) 
CALL CLOSE (I DAT I ) 
CALL KDPOP 
RETURN , 


9999 CONTINUE 

CALL CLOSE (I DAT I) 
CALL CLOSE (IDATJ) 


ERROR RETURN 


194 








CTXINPT T-X-I-N-P-T 

C 


c 

C ASCII I/O FOR THE TEXTURE PROGRAMS 


C 

c 

c 

c 

c 

c 

c, 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

i; 


PROGRAM TITLE 

VERSION 

AUTHOR 

DATE 

UPDATE 

PROGRAM LANGUAGE 
IMPLEMENTED ON 
DOCUMENTED BY 
PURPOSE 


TXINPT 

A 

CH IN-HUANG CHEN 
FEBRUARY 1975 

FORTRAN IV 
PDF' 15 

CH IN-HUANG CHEN 


THIS ROUTINE GETS THE NECESSARY PARAMENTERS FOR THE 
TEXTURE TRANSFORM PACKAGE 

ENTRY POINT TXINPTCNFUNC, NDIS. FILNMP, FILNMQ, IBOUT. 


PCLCT ) 


ARGUMENT LIST 

NFUNC PARAMETER USED TO DETERMINE WHICH FUNCTION 

COMPUTES THE JDATA IMAGE 
NFUNC=1 FOR SUM PROBABILITY FEATURE 
NFUNC=2 FOR ANGULAR MOMENTUM FEATURE 

NFLINC=3 FOR ENTROPY FEATURE 

NFUNC=4 FOR GRADIENT FEATURE 

NFUNC=5 FOR NORMALIZED ARRAY WHICH HAS 

BEEN EQUAL PROBABILITY QUANTIZED 
NDIS SPATIAL DISTANCE TO BE USED TO GENERATE LEX ARRAYS 

FILNMP INPUT FILE NAME 

FILNMQ OUTPUT FILE NAME 

I BOUT ERROR MESSAGE OUTPUT . DAT SLOT 

PCLCT PERCENT OF LINES COUNTED IN GENERATING THE 

FOUR NEIGHBOR GRAY TONE MATRICES (LEX ARRAYS) 


SUBROUTINE TX INPT < NFUNC, NDIS, FILNMP, FILNMQ, lEOUT, PCLCT) 
DOUBLE INTEGER FILNMP, FILNMQ, FDATE 
DIMENSION FILNMP(2) ,FILNMQ<2), FDATE(S) 

TOUT = 6 
IDIN =4 

GET PARAMETERS 


200 WRITECIOUT, 100) 

100 FORMAT (••■ TYPE NFUNC, NDIS, I BOUT, PCLCT, I/O FILE NAMES’) 
WRITE( lOUT, 110) 

110 FORMAT(- (FORMAT IS 315, F4. 2, A'9, A9 ) - ) 

READdDIN, 101 ) NFUNC, NDIS, IBOUT, PCLCT, FILNMP, FILNMQ 

101 F0RMAT(3I5, F4. 2, A5, A4, A5, A4) 

WR ITE ( lOUT, 1 02 ) NFUNC, ND I S, I BOUT, PCLCT, F I LNMP, F I LNMQ 
IF( IBOUT. NE. lOLlT) WRITE/ IBOUT, 102)NFUNC, NDIS, IBOUT, FILNMP, 
2FILNMQ 

102 FORMAT( IX, 315, 2X, F4, 2, 2X, A5, A4, 2X, AS, A4) 

C ■ 

CALL ADATE( FDATE) 

WRITE( TOUT, 405) FDATE 


196 


IF( IBOUT. NE. IDUT)WRITE( IBOUT, 405) FDATE 
i 405 FORMAK IX, 3A5) 

I- RETURN 

END 

! 

i 

I 
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CTX.JDM 


T-X-J-D-M 


C 

C TEXTURE JDATA MAINLINE 

C 


TXJDM 

A 

ROBERT M HARALICK 
NOVEMBER 1974 
FEBRUARY 1975 
CHIN-HUANG CHEN 
FORTRAN IV 
POP 15 


C 

C THIS ROUTINE IS THE MAIN LINE FOR THE JDATA GENERATION 

C DISPLAY . INPUT PARAMETERS ARE FUNCTION TYPE, SPATIAL DISTANCE 

C RELATIONSHIP, ERROR MESSAGE OUTPUT .DAT SLOT, PERCENT OF LINES 

C COUNTED IN GENERATING THE FOUR NEIGHBOR GRAY TONE MATRICES, 

C INPUT FILE NAME, AND OUTPUT FILE NAME. 

C 

C SUBROUTINES CALLED: 

C TXINPT 

C ERROR 

C TXTMN 

C SDKINL 

C • SKPDSC 

C FPLXIT 

C SREAD 

C INDEX 

C FUNCl 

C FUMC2 

C FUNC3 

C FUNC4 

C FUNC5 

C EQPQNT 

C LEXEQP 

C PLXIT 

C INDEX 

C SREAD 


C PROGRAM TITLE: 

C VERSION: 

C AUTHOR: 

C DATE 

C UPDATE 

C 

C PROGRAM LANGUAGE: 

C IMPLEMENTED ON: 

C PURPOSE: 


C 

DOUBLE INTEGER NPPCAL,NTOTAL, FI LNMP, FI LNMQ 
DIMENSION FILMMGK 2) , FXLNMPX 2) 

COMMON IWORKX 8000 ), IWRKC 7000) 

COMMON ./DFA./‘ NG, F(50) 

COMMON./DFB./AMEAN, VAR, NPPCAL, NTOTAL, START, END, NCALL, NNTERS, DANGE 
COMMON 10/ N:SIZ E, I DUM( 2048) 

COMMON ./TXT/ ITXT< 10) 

DATA lOUT, IDATN, IDATQ, NDIM/6, 2, 1, 15000/ 

NSIZE = 1024 
200 CONTINUE 

C 

C GET INPUT PARAMETERS 


198 


CAll. TXIMPTCMFUlvIC, NDIS, FILMMP, FILNIMQ, IBOUT, PCLCT) 


WRITE THE INPUT IMAGE IDENTIFICATION 
BLOCK 

CALL LST I D ( I D ATK> F I LNMP , I BOUT ..1,1 EV , @304 ) 

COMPUTE THE INTEGER TEXTURE IMAGE 

CALL TXTMNC I WORK, NDIM, FILNMP, FILNMQ, NDIS, NFUNC, 

2NXMIN, MXMAX, PCLCT, lEV, @310) ' 

WRITE THE OUTPUT IMAGE IDENTIFICATION 
BLOCK 

CALL LSTID<IDATQ, FILNMQ, IB0UT,1, lEV, ©304) 


CALL CLOSE (I BOUT) 

GO TO 200 
04 IERR-1 

GO TO '500 
10 IERR=2 

00 CALL ERROR < I ERR, lEV, I OUT, I BOUT) 

GO TO 200 


PROGRAM 

TITLE 

PROGRAMMER 


UPDATE 

DOCUMEM- 
TAT ION 

eOMPUTER 

REQUIRED 

PROGRAM 

LANGUAGE 

PURPOSE 

METHOD 

ENTRY POINT 
ARGUMENTS 


SUBROUTINE TXTMN 


MODIFIED BY A. SINGH, OCTOBER hl'l 

MODIFIED FOR PDF' BY ROBERT M HARALICK 5/1/73 

7/10/73 

2/2/74 

6/30/74 

GE MONAGHAN 9/20/74 

RM HARALICK 10/10/74 

CH IN-HUANG CHEN 2/22/75 

A, SINGH, OCTOBER S<72 


PDF'- 15 


FORTRAN IV 


TXTMN IS THE MAINLINE SUBROUTINE FOR THE TEXTURE 
ROUTINE PACKAGE TO COMPUTE THE JDATA IMAGE. 

TXTMN DOES THE FOLLOWING - 

TAKES IN LABELS AND PARAMETERS FROM ARGUMENT LIST, 
READS IN THE IMAGE FROM FILE (02), 

SETS THE MAXIMUM AND MINIMUM GRAY LEVELS, 
SETS UP DYNAMIC ALLOCATION OF PARAMETERS AND 
CALLS THE REST OF THE SUBROUTINES. 

CALL TXTMN( IWORK, NDIM, S, T, NDIS, IFUNC, 

NXMIN, NXMAX, PCLCT, lEV, ERRl) 

SCRATCH ARRAY WHERE THE IMAGE IS READ IN 
AND THEN LATER IT IS USED FOR DYNAMIC 
ALLOCATION. 

SIZE OF SCRATCH ARRAY, NDIM SHOULD BE 
EITHER NUMPPL-is-NUMLIN OR 

2-«-(M-K-(A+l)+l)+4-«-B-«-<B+l), WHICH EVER ONE 
IS LARGER. 

A IS THE NUMBER OF PO I NTS /LINE IN THE 
IMAGE, B IS THE MAXIMUM NUMBER OF GRAY 
LEVELS POSSIBLE AND M IS THE LARGEST 
REDUCTION DISTANCE THE PROGRAM WILL RUN 
WITH. 

NAME OF FILE THE IMAGE IS ON 

NAME OF FILE WHERE THE JDATA IMAGE IS CREATED 
SPACING BETWEEN NEIGHBORLY RESOLUTION CELLS 
PARAMETER USED TO DETERMINE WHICH FUNCTION 
COMPUTES THE JDATA IMAGE. 

IFUNC=1 FOR SUM PROBABILITY 
IFUNC=2 FOR ANGULAR MOMENTUM FEATURE 


IWORK 


NDIM 


T 

NDIS 

IFUNC 


I 

I c 

I 

i c 

] 

I c 

J 

I c 


I 

i C 

: c 


PARAMETERS 
AND ARRAYS 


IFUNC=3 FOR ENTROPY FEATURE 
IFUMr;"4 FOR GRADIENT FEATURE OF THE IMAGE 
IFlJNr '-0 FOR NORMALIZED LEX ARRAY WHICH HAS BEEN 
EQUAL PROBABILITY QUANTIZED 
NXMIN IS THE MINIMUM ON THE JDATA IMAGE. 

MX MAX IS THE MAXIMUM ON THE JDATA IMAGE 

PCLCT IS THE PERCENT OF LINES COUNTED 
lEY INTEGER EVENT VARIABLE 

ERR I ALTERNATE ERROR RETURN 


NUMLIN NUMBER OF LINES IN THE INPUT IMAGE 

NUMPPL NUMBER OF POINTS PER LINE IN THE INPUT IMAGE 

I MAX MAXIMUM GRAY LEVEL 

I MIN MINIMUM GRAY LEVEL 

LEAST! =IMIN-1 

NOBL NUMBER OF GRAY LEVELS 

NBUBL =N0BL^^-(N0BL+1 )/2 IS THE SIZE OF A LEX ARRAY 

NIDATA, NJDATA, NLEXl. 2, 3, 4, NFUNC AND NTOT ARE POINTERS 
FOR DYNAMIC ALLOCATION IN I WORK. 


INPUT AND OUTPUT 


! C 
' c 
i c 
i c 


READ IN FROM FILE (02) 

ERROR FOR INCORRECT SIZE OF I WORK, 
ERROR IF PARAMETER IFUNC HAS BEEN 
INITIALIZED INCORRECTLY. 

INPUT IMAGE ON FILE CODE IDATI. 

WRITES JDATA IMAGE ON FILE CODE IDATJ. 


SUBPROGRAMS 

REQUIRED 

RETURNS 


FPLXIT, PLXIT, INDEX.. FUNCl, F1JNC2, FUNC3, FUNC4, FUNC5 
SEEK SYSTEM LIBRARY 

CLOSE SYSTEM LIBRARY 


NORMAL AND ALTERNATE 
PROGRAM TERMINATED FOR INCORRECT 
S I ZE OF I WORK, ERROR I F I FUNC 
INITIALISED INCORRECTLY. 


CALLED BY 


MAIN LINE PROGRAM TXJDM 


SUBROUTINE TXTMN( IWORK, NDIM, S, T, NDIS, IFUNC, NXMIN, 
2NXMAX, PCLCT, lEV, ERR! ) 

INTEGER ERRl 
DOUBLE INTEGER FUND 
DOUBLE INTEGER FDATE, S, T, A, B 
DOUBLE INTEGER L.EX, Cl, R1, CC!,RR1 




- * p ■* 


. ■ 


DIMENSION I WORK (1 ), LEXL( 1), LEX2( 1) , LEX3U) , LEX4( 1>, IDATA( 1, 1), T(2) 
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D I MENS I ON I DENT < 20 ) , S ( 2 ) ; JDENT < 20 ) , FDATE < 3 ) 

D I MENS I ON JDAT A ( 1 , 1 ) , I PT ( 1 ) , FUNC < 1 , 1) 

DIMENSION CCKS), RRl (8) , LEX < 13 ) , OK 2 ) , Rl(2) 


I DATA < NUMPPL, MM 1 ) , JDATA < NLIMPPL. MMl ) , LEX 1 ( NBUEL ) , LEX2 < NBUBL ) . 
LEX3 ( NBUBL ) , LEX4 ( NBUBL ) , IPT ( MMAX ) , FUNOT < NBUBL. 4 ) 


COMMON /TXT/ IMAX. IMIN, NUMPPL, NUMLIN. NBUBL. NOBL. LEASTl 
DATA A, B, IZ. lONE, ITWO/--TXTMN-‘ , - - , 0, 1, 2/ 

DATA I DAT I. I DAT J/2, 3/ 


DATA LEX / LEX - , - ARRAY . 1 1 
DATA Cl, R1/-C0L •■. •ROW ■■ 

DATA LC1/ ‘ L- 1 " . ' U2’ . • Lc'"^ , " C4‘ 
DATA RR 1 / R 1 . - R2 R3 -■ . R4 


■■■ / 

. ■ L5 ' . ' C6>‘ . ' L"/ ' . 
, •■•R5'. •R6'. -'R?-'. 


CS •' / 
•••RS-/ 


CALL KDPUSH(A. B') 

CALL SDKINL< IDATI. S. IDENT, 1, lEV. ERRl ) 
NUMPPL=IDENT(6) 

NUMLIN=IDENT(7) 

IMIN=IDENT( 15) 

IMAX=IDENT( 16) 


LEAST1=IMIN-1 
NOBL= I MAX-LEAST 1 
NBUBL=NOBL-«^ < NOBL+1 ) /2 


SET DYNAMIC ALLOCATION PARAMETERS 

SINCE THE SIZE OF IDATA, JDATA AND IPT ARE DIFFERENT FOR DIFFERENT 
REDUCTIONS, THE MAXIMUM SPACE THEY WILL REQUIRE HAS TO BE 
RESERVED. IDATA AND JDATA THIS WILL BE < NDIS+1 ) *NUMPPL, AND FOR 
IPT JUST NDIS+1. 

MM1=NDIS+1 

NIDATA=1 

N JDAT A=N I D AT A+NUMPPL-tt-MM 1 
NLEX 1 =N JD ATA+NUMPPLi^MM 1 
NLEX2=NLEX 1+NBUBL 
NLEX3=NLEX2+NBUBL 
NLEX4=NLEX3+NBUBL 
NFUNC=NLEX1 
N I PT= ( NLEX 4+NBUBL ) *2 
NT0T=NIPT+MM1 

CHECKING IF THE SIZE OF IWORK IS ENOUGH 

IF(NTOT. GT. NDIM)GO TO 78 
NBIG=NDIM-NTOT T 
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Ar:i. ,II JST THE D I HENS I QMS 

C A1 ...I .. AD J2 ( I DATA , I WORK ( N I DATA ) , MUMPPL ) 

CALL ADJ2 < .JDATA, I WORK < NJDATA ) , NUMPPL ) 

CAL L AD J J. ( LEX 1. , I WORK < NLEX 1 ) ) 

CALL ADJULEX2; IW0RK<WLEX2) ) 

CALL ADJHL.EX3. IWORK(NLEXS) ) 

CAL.L ADJ 1 ( LEX 4. I WORK ( NLEX4 ) ) 

CAI-„L AD J2 < FUNC. I WORK •: MFUL4C ) , NBUBL ) 

CALL AD JiaPT.IWORK< NIPT)) 


ZERO OUT THE SCRATCH AREA 


DO 30 JLK=1. NDIM 
0 IW0RK(JLK)=0 

SKIP THE DESCRIPTOR RECORDS 
CALL SKPDSC ( I DAT I , I DENT. I EV, ERR 1 ) 


COMPUTE THE FOUR LEX ARRAYS 

CALL FPLXrn IDATI. TDATA, LEX 1 , LEX2, LEX3. LEX4.. IPT, IDENT, MMl, 
2PCLCT, lEU, ERRl ) ' 


WRITE OUT THE LEX ARRAYS 


CALL 

r:At...L 

CALI, 

CAL.L 


IMTRXFULEXI, 8, 8. 8. 
IMTRXPCLEX2/ 8. 8. 8. 


IMTRXPXLEXS, 8> 
IMTRXP( LEX4, 8, 



LEX, Cl/Rl; CCl, RRL 4) 
LEX, Cl, Rl, CCl, RRl, 4) 
LEX, Cl, Rl, CCl, RRl, 4) 
LEX, Cl, Rl, CCl, RRl, 4) 


CALL PROPER FUNCTION; SUBPROORAM 

IEV--5011 

IF( IFUWC. ECl T ) CALL FUNCl (LEXl, LEX2, LEX3, LEX4, FUNG, NBUBL) 

IF( IFUNC. EQ. 2) CALL FUWC2XLEX1, LEX2, LEX3, LEX4, FUNC, NBUBL) 

TF( IFUNC. EG! 3) CALL FUWC3':LEX1, LEX2, LEX3, LEX4, FUNC, NBUBL) 

X F ( I FUNG. EQ. 4 ) CALL FUWC4 (LEXl , LEX2, LEX 3, LEX 4 , FUNC, NBUBL ) 

I.F( IFUNC. EQ. 5) CALL FUNC5( LEX 1, LEX2, LEX3, LEX4, FUNC, NBUBL) 

IF( IFUNC. LE. 0. OR. IFUNC. GE. 6) RETURN ERRl 
DO 78 1^1, 20 
JDENT( r)=rDENT( I) 

CONTINUE 


JDENT(5)=10 
JDEWT( 18) “1 
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JDENT(10)=3 
JDENT(11)=512 
JDENT< 15)=-256 
JDENK 16)=255 

CALL C:PYDSC< IDATI> S; IDATJ, T. JDENT, lEV, ERRl ) 

CALL ADATE(FDATE) 

NW=JDENT( 12)*2 

WRITE< IDATJ) A, E:. FDATE, S. ( IZi I = 15> NW) 

WRITE(IDATJ) lONE, < 12, 1=2, NW) 

WRITE< IDATJ) ITWO, IFUNC, MDIS, < IZ, 1=3, NW) 

CALL PLXIK IDATI, IDATJ, IDATA, JDATA, IDENT, FUNC, IPT, NBUBL, MMl, 
NXMIN, NXMAX, JDENT, lEV, ERRl ) 

CONTINUE 
CALL KDPOP 
RETURN 

ERROR RETURN FOR NOT ENOUGH WORK SPACE 

IEV=-5010 
RETURN ERRl 


